AD-A267  666 


SUBJECTED  TO 

THERMOMECHANICAL  CYCLIC  LOADING 


Demirkan  Coker  and  Noel  E.  Ashbaugh 


University  of  Dayton  Research  Institute 
300  College  Park 
Dayton,  OH  45469-0001 


December  1992 

Final  Report  for  01/01/92-12/01/92 


‘Tr1 1  .*■> 

fi  V  3  8 

«  s',  r  V  4 


^AUGO  4 1993,1  § 


Approved  for  Public  Release;  Distribution  is  Unlimited 


Materials  Directorate 
Wright  Laboratory 
Air  Force  Materiel  Command 
Wright  Patterson  AFB  OH  45433-7734 

'  93-17420 

vV  | 

!  . .  '  '■*  vO  a 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used  for 
any  purpose  other  than  in  connection  with  a  definitely  Government-related 
procurement,  the  United  States  Government  incurs  no  responsibility  or  any 
obligation  whatsoever.  The  fact  that  the  government  may  have  formulated  or 
in  any  way  supplied  the  said  drawings,  specifications,  or  other  data,  is  not 
to  be  regarded  by  implication,  or  otherwise  in  any  manner  construed,  as 
licensing  the  holder,  or  any  other  person  or  corporation;  or  as  conveying 
any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention 
that  may  in  any  way  be  related  thereto. 


This  report  is  releasable  to  the  National  Technical  Information  Service 
(NTIS).  At  NTIS,  it  will  be  available  to  the  general  public,  including 
foreign  nations. 


This  technical  report  has  been 

tion. 


Project  Engineer 
Materials  Behavior  Branch 


reviewed  and  is  approved  for  publica- 


ALLAN  W.  GUNDERSON,  Chief 
Materials  Behavior  Branch 
Metals  and  Ceramics  Division 


NORMAN  M.  GEYER 
Acting  Deputy  Chief 
Metals  and  Ceramics  Division 
Materials  Directorate 


If  your  address  has  changed,  if  you  wish  to  be  removed  from  our  mailing 
list,  or  if  the  addressee  is  no  longer  employed  by  your  organization  please 
notify  WL/MLLN  .  WPAFB,  OH  45433-  7817  to  help  us  maintain  a  current 
mailing  list. 


Copies  of  this  report  should  not  be  returned  unless  return  is  required  by 
security  considerations,  contractual  obligations,  or  notice  on  a  specific 
document. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704  0188 


Public  reportinq  burden  for  thu  collection  of  information  is  estimated  to  average  t  Hour  per  response  including  the  time  tor  reviewing  instructions,  searchmq  easting  data  sources, 
gathering  and  maintaining  the  data  needed,  and  completing  ana  reviewing  the  collection  of  information  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 
collection  of  information,  including  suggestions  for  reducing  this  Durden  to  Washington  Headauarters  Services.  Directorate  for  information  Operations  and  Reports,  1215  Jefferson 
Davis  Highway.  Suite  1204.  Arlington,  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  (0  704-0186).  Washington,  DC  20503 


AGENCY  USE  ONLY  (Leave  blank) 


2  REPORT  DATE 

December  1992 


REPORT  TYPE  AND  OATES  COVERED 

Interim:  January  1992-December  1992 


4.  TITLE  AND  SUBTITLE 

Elastic-Plastic  Finite-Difference  Analysis  of 
Unidirectional  Composites  Subjected  to  Thermomechanical 
Cyclic  Loading 


6.  AUTHOR(S) 

Demirkan  Coker  and  Noel  E.  Ashbaugh 


5.  FUNDING  NUMBERS 

F33615-91-C-5606 

PE-61102 

PR-2302 

TA-P1 

WU-03 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Dayton  Research  Institute 
300  College  Park 
Dayton,  OH  45469-0001 


B.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Materials  Directorate 

Wright  Laboratory 

Air  Force  Materiel  Command 

Wright-Patterson  Air  Force  Base,  OH  45433-7734 


10.  SPONSORING  /  MONITORING 
AGENCY  REPORT  NUMBER 


WL-TR-93-4043 


12a.  DISTRIBUTION  /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  public  release;  distribution  is  unlimited. 


13.  ABSTRACT  (Maximum  ?00  words) 

An  analytical  tool  was  developed  to  model  a  unidirectional  composite  subjected  to 
thermomechanical  cyclic  loading  and  processing  conditions.  The  finite  difference  method  was 
incorporated  into  a  PC  compatible  computer  code,  FIDEP  (Finite-Difference  Code  for  Elastic-Plastic 
Analysis  of  Composites).  FIDEP  provides  an  efficient  numerical  procedure  for  analyzing  a  variety  of 
problems  involving  thermal  and  mechanical  cycling.  The  procedure  allows  the  modeling  of  the 
constituent  materials  as  elastic-plastic  with  temperature  dependent  properties.  The  concentric  cylinder 
approximation  allows  the  computations  to  capture  the  three-dimensional  aspects  of  the  stress  state  in 
a  real  composite.  Results  for  a  thermal  cool-down  in  an  SCS-6/Ti-24AI-1 1  Nb  unidirectional  composite 
compare  well  with  those  obtained  using  the  finite  element  method.  Several  problems  were  solved  for 
thermomechanical  loading  conditions  and  demonstrate  the  three-dimensional  nature  of  the  stress  fields 
in  the  matrix  material. 


14.  SUBJECT  TERMS 


Concentric  cylinder  model,  e'astic  plastic,  finite  difference  methods,  metal 
matrix  composite,  thermomechanical  fatigue,  unidirectional  composite. 


1  SECURITY  CLASSIFICATION 
OF  REPORT 

Unclassified 


I.  SECURITY  CLASSIFICATION 
OF  THIS  PAGE 

Unclassified 


19.  SECURITY  CLASSIFICATION 
OF  ABSTRACT 

Unclassified 


16  PRICE  CODE 


20.  LIMITATION  OF  ABSTRACT 


NSN  7540  01 -.280-5500 


Standard  Form  298  (Rev  2*89) 

present**  by  ANSI  Std  239 
298-102 


TABLE  OF  CONTENTS 


LIST  OF  FIGURES . 

LIST  OF  TABLES . 

FOREWORD . 

1.  INTRODUCTION . 

2.  BACKGROUND . 

3.  THEORY  . 

3.1  Concentric  Cylinder  Model . 

3.2  Governing  Equations  . 

3.3  Computation  of  Axial  Strain  . 

3.4  Finite  Difference  Formulation  . 

3.5  Plasticity  Theory . 

3.5.1  Prandtl-Reuss  Relations  . 

3.5.2  Plastic  Strain-Total  Strain  Plasticity  Equations 

3.6  Solution  Technique . 

4.  IMPLEMENTATION  IN  FIDEP . 

4.1  Algorithm  . 

4.2  Faster  Convergence  Scheme . 

5.  PROGRAM . 

5.1  General  Description  . 

5.2  Program  Operation  . 

5.2.1  Input  . 

5.2.2  Output . 


5.2.3  Program . 32 

5.2.4  Execution  of  FIDEP  on  a  VAX/VMS  Machine . 32 

6.  DEMONSTRATION  PROBLEMS . 34 

6. 1  Material  Properties  . 34 

6.2  Comparison  with  Elastic  Solution . 34 

6.3  Comparison  with  Finite  Element  Method  . 38 

6.3.1  Cool-Down  from  Processing  Temperature . 38 

6.3.2  Cyclic  Loading . 45 

6.4  Thermo-Mechanical  Fatigue  (TMF) . 48 

REFERENCES . 57 

APPENDIX  1.  LISTING  OF  FIDEP  SOURCE  CODE . 60 

APPENDIX  2.  ELASTIC  SOLUTION  OF  TWO  CONCENTRIC 

CYLINDERS . 84 

APPENDIX  3.  SAMPLE  OUTPUT  FILES . 85 


IV 


LIST  OF  FIGURES 


Figure  3. 1  Concentric  Cylinder  Idealization  of  a  Unidirectional  Composite  and 

the  Concentric  Cylinder  Model . 5 

Figure  3.2  Discretization  of  the  Concentric  Cylinder  Model . 10 

Figure  3.3  The  Linear  System  of  Equations  Resulting  from  the  Finite 

Differences  Formulation  of  the  Governing  Equations  for  the 
Concentric  Cylinder  Model . 14 

Figure  3.4  Schematic  of  the  Iterative  Technique  for  Computing  Plastic  Strains 

Using  the  Prandtl-Reuss  Relations . 19 

Figure  3.5  Schematic  of  the  Iterative  Technique  for  Computing  Plastic  Strains 

Using  the  Modified  Prandtl-Reuss  Relations . 19 

Figure  4. 1  Elastic-Plastic  Algorithm  Used  in  FIDEP . 22 

Figure  4.2  Pseudo-Code  for  the  Computer  Program  FIDEP . 23 

Figure  4.3  Schematic  for  the  Computation  of  a  New  Estimate  for  the  Plastic 

Strains  Using  the  Previous  Estimates  and  their  Differences . 26 

Figure  5.1  Example  Loading  History  Input  File,  FDLOAD.DAT . 29 

Figure  5.2  Example  Material  Properties  Input  File,  FDMAT.DAT . 29 

Figure  5.3  Example  Command  File  for  Running  Batch  Jobs  on  a  VAX/VMS 

Computer  System . 33 

Figure  6. 1  Variation  of  Elastic  Stresses  with  Radius  Obtained  by  Elastic 

Closed-Form  Solution  and  FIDEP  for  AT=-100°C . 36 

Figure  6.2  Elastic  and  Elastic-Plastic  Predictions  of  Stresses  at  a  Cross-Section 
of  the  Composite  after  Cool-Down  from  the  Processing  Temperature 
of  1010°C . 37 

Figure  6.3  a)  Loading  Input  File,  and  b)  Schematic  of  the  Temperature  History 
for  a  Thermal  Cool-Down  Simulation  of  SCS-6/Ti-24Al- 1 1  Nb 
Unidirectional  Composite . 39 

Figure  6.4  Finite  Element  Model  and  the  Boundary  Conditions  for  an 

Axisymmetric  Concentric  Cylinder  Geometry  Under  Generalized 

Plane  Strain  Condition . 40 

Figure  6.5  Finite  Element  and  FIDEP  Predictions  for  the  Stresses  at  the 

Interface  in  Ti-24A1-1  INb  Matrix  after  Cool-Down  from  Processing 
Temperature  of  1010°C . 41 


Figure  6.6  Finite  Element  and  FIDEP  Predictions  for  the  Plastic  Strains  at  the 
Interface  in  Ti-24A1-1  INb  Matrix  after  Cool-Down  from  Processing 


Temperature  of  1010°C . 42 

Figure  6.7  Variation  of  the  Stresses  Across  the  Cross-Section  of  the  CCM 

After  Cool-Down  from  Processing  Temperature  of  101  CPC . 43 

Figure  6.8  Variation  of  the  Plastic  Strains  Across  the  Cross-Section  of  the 

CCM  After  Cool-Down  from  Processing  Temperature  of  1010°C . 44 

Figure  6.9  a)  Loading  Input  File,  and  b)  Schematic  of  the  Temperature  History 
for  Thermal  Cyclic  Loading  Simulation  of  SCS-6/Ti-24Al-l  INb 
Undirectional  Composite . 46 

Figure  6. 10  FIDEP  and  FEM  Predictions  for  the  Effective  and  Radial  Stresses  in 
the  Matrix  at  the  Fiber-Matrix  Interface  for  Thermal  Loading 
History  shown  in  Fig.  6.9 . 47 

Figure  6. 1 1  a)  Loading  Input  Files,  and  b)  Schematic  of  the  Temperature  and 

Stress  History  for  Typical  In-Phase  and  Out-of-Phase  TMF  Cycles, . 49 

Figure  6.12a  Variation  of  Predicted  Stresses  with  Radius  for  an  Out-of-Phase 

TMF  Cycle  at  150°C . 50 

Figure  6. 12b  Variation  of  Predicted  Stresses  with  Radius  for  an  Out-of-Phase 

TMF  Cycle  at  650°C . 51 

Figure  6.13  Axial  Stress  Predictions  in  the  Ti-24A1-1  INb  Matrix  at  the  Fiber- 
Matrix  Interface  for  In-Phase  and  Out-of-Phase  TMF  Cyclic 
Loading . 53 

Figure  6. 1 4  Axial  Stress  Peaks  Predicted  in  the  SCS-6  Fiber  for  In-Phase  and 

Out-of-Phase  TMF  Cyclic  Loading . 54 

Figure  6.15  Radial  and  Hoop  Stress  Predictions  in  Ti-24A1-1  INb  Matrix  at  the 

Fiber-Matrix  Interface  in  TMF  Cyclic  Loading . 55 


v; 


LIST  OF  TABLES 


Table  5.1  Format  for  the  Input  Loading  File,  FDLOAD.DAT . 30 

Table  5.2  Format  for  the  Material  Properties  File,  FDMAT.DAT . 30 

Table  6. 1  Mechanical  Properties  for  SCS-6  Fiber  and  Ti-24A1- 1 1  Nb  Matrix . 35 


FOREWORD 


This  report  documents  a  computer  program  that  was  developed  as  part  of  an  investigation  of 
the  mechanical  behavior  of  metal  matrix  composites.  The  investigators  were  Demirkan 
Coker  and  Noel  E.  Ashbaugh  of  the  Structural  Integrity  Division,  University  of  Dayton 
Research  Institute,  Dayton  OH.  The  research  was  conducted  at  the  Materials  Behavior 
Branch,  Metals  and  Ceramics  Division,  Materials  Directorate,  Wright  Laboratory 
(WL/MLLN)  Wright-Patterson  AFB  OH,  under  Contract  No.  F33615-91-C-5606.  The 
contract  was  administered  under  the  direction  of  WL/MLLN  by  Mr.  Jay  R.  Jira.  Dr.  Noel 
E.  Ashbaugh  was  the  Principal  Investigator  and  Dr.  Joseph  P.  Gallagher  was  the  Program 
Manager. 


CHAPTER  1 
INTRODUCTION 


Interest  in  metal  matrix  composites  (MMC)  for  high-temperature  aerospace 
applications  has  grown  in  recent  years.  Because  of  the  high  use  temperatures  envisioned, 
proposed  applications  of  MMCs  almost  always  involve  both  thermal  and  mechanical  cyclic 
loading.  Consequently,  the  accurate  prediction  of  the  thermomechanical  fatigue  (TMF)  life 
of  these  components  is  a  critical  aspect  of  the  design  process.  Because  of  the  mismatch  in 
coefficient  of  thermal  expansion  between  fiber  and  matrix,  internal  thermal  stresses  are 
produced  during  cool-down  from  the  processing  temperature  and  subsequent  thermal 
cycling.  These  stresses  must  be  accounted  for  in  the  analysis  in  addition  to  those  produced 
by  mechanical  loading.  Prediction  of  the  mechanical  behavior  and  life  of  a  composite 
material,  therefore,  requires  a  knowledge  of  the  individual  components  of  stress  or  strain  in 
the  fiber  and/or  matrix  as  well  as  the  overall  applied  stresses. 

In  this  investigation,  an  axisymmetric  concentric  cylinder  geometry  was  used  to 
model  the  unidirectional  composite  in  which  the  fiber  and  matrix  were  represented  by  the 
core  cylinder  and  outer  cylinder,  respectively.  The  concentric  cylinder  geometry  allowed  for 
a  relatively  simplified  analysis  and  was  adequate  for  predicting  the  complex  three- 
dimensional  stress  distributions  in  the  fiber  and  the  matrix.  The  analysis  allowed  for  axial 
and  radial  loading  and  more  than  two  constituents  with  bilinear  elastic-plastic  properties. 
The  analysis  also  accounted  for  thermomechanical  loading  and  temperature  dependent 
material  properties.  These  aspects  were  incorporated  into  a  FORTRAN  program,  FIDEP 
(Finite  Difference  Code  for  Elastic-Plastic  Analysis). 

This  report  summarizes  the  derivation  and  application  of  the  concentric  cylinder 
model  to  predict  the  micromechanical  stresses  in  an  SCS-6/Ti-24Al-l  INb  metal  matrix 
composite.  The  next  chapter  summarizes  some  of  the  past  work  dune  on  micromechanical 
modeling  of  composites.  Chapter  3  describes  the  theory  for  the  concentric  cylinder  model 
and  summarizes  plasticity  theory.  In  Chapter  4  the  implementation  of  the  theory  to  the 
computer  program,  FIDEP,  is  described  and  the  procedures  for  using  this  program  are 
explained  in  Chapter  5.  Finally  some  example  runs  are  conducted  with  this  program  and 
the  results  are  presented  . 


CHAPTER  2 
BACKGROUND 

Various  appi caches  have  been  adopted  in  the  past  to  calculate  the  micromechanical 
stresses  in  the  fiber  and  the  matrix  due  to  both  thermal  and  mechanical  loading.  The  most 
widespread  approach  involves  use  of  the  finite  element  method  to  model  a  representative 
vo’  .me  element  (RVE)  and  compute  the  stress  distributions  around  the  fiber.  Common 
RVE  geometries  that  have  been  used  include  either  a  square  or  a  cylindrical  fiber  in  a  square 
matrix  and  a  cylindrical  fiber  in  a  cylindrical  matrix.  Different  constitutive  behaviors 
include  elastic-plastic  behavior  [1-3],  time-dependent  behavior  using  a  classical  creep  law 
[4]  and  unified  viscoplastic  theories  such  as  the  Bodner-Partom  model  [5],  and  Bodner- 
Partom  with  backstress  [6].  However,  finite  element  codes  are  sually  run  on  mainframe 
computers,  are  time  consuming,  allow  for  solution  of  only  one  problem  at  a  time,  and 
require  familiarity  with  finite  element  analysis  and  the  code.  To  conduct  a  parametric  study 
and  to  support  the  design  of  a  large  number  of  experiments  on  different  composite  systems, 
a  more  practical  approach  is  desirable. 

One  simple  approach  to  micromechanical  modeling  of  composites  has  been  to  use 
average  stresses  in  the  constituents.  This  method  has  the  advantage  of  being  easy  to 
implement  into  a  code  and  of  easily  being  extended  to  analyze  laminates.  Approaches  using 
this  method  include  analysis  of  square  fiber  in  a  square  matrix  and  the  vanishing  fiber 
diameter  (VFD)  model  developed  by  Dvorak  and  Bahei-El-Din  [7].  In  the  VFD  model,  the 
presence  of  the  fibers  is  assumed  not  to  influence  the  transverse  stresses.  The  transverse 
properties  of  a  lamina  are  easily  computed  using  this  assumption  while  the  longitudinal 
properties  are  calculated  from  the  rule  of  mixtures.  This  analysis  can  be  carried  out  for  an 
elastic  fiber  surrounded  by  a  matrix  material  having  a  variety  of  constitutive  relations  to 
calculate  ply  properties  which,  in  turn,  can  be  incorporated  into  lamination  theory  to  analyze 
a  composite.  Such  a  procedure  has  been  implemented  into  computer  codes  such  as 
AGLPLY  1 8 1  which  uses  thermoplastic  matrix  behavior  and  V1SCOPLY  |9|  which 
incorporates  the  Eisenberg  creep  model  to  predict  average  stresses  in  the  constituents  and  in 
a  symmetric  composite  laminate  subjected  to  thermomechanical  loading.  Aboudi  ]  10] 
modeled  a  square  fiber  in  a  square  matrix  subcell  using  first  order  displacement  expansions. 
The  unified  theory  of  Bodner  and  Partotn  was  used  to  compute  inelastic  strains.  Hopkins 
and  Chamis  ( 1 1 1  and  Sun  1 12]  conducted  strength-of-materials  type  analyses  for  a  square 
fiber  in  a  square  matrix  cell  to  obtain  expressions  for  the  constituent  microstresses.  Chamis 
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and  Hopkins  [13]  modified  these  expressions  based  on  exneiimental  data  for  uniaxial 
lamina  and  three-dimensional  computer  simulations  of  composite  behavior.  The  results 
were  incorporated  into  a  computer  program,  METCAN  [14],  which  treats  material 
nonlinearity  at  the  constituent  level  where  material  behavior  is  modeled  using  a  time- 
temperature-stress  dependence  of  constituent  properties. 

The  main  drawback  to  all  of  these  simplified  material  models  is  the  assumption  of 
average  uniform  stress  in  the  representative  volume  element.  This  approach  fails  to  take  into 
account  the  triaxial  stresses  that  arise  due  to  the  mismatch  of  the  Poisson's  ratio  and  the 
coefficient  of  thermal  expansion  of  the  constituents.  A  detailed  review  of  these  codes  and 
comparisons  of  the  stresses  using  these  methods  and  using  finite  element  analysis  is  given 
by  Bigelow  et  al.  [15]. 

For  a  more  accurate  analytical  treatment  of  the  triaxial  stresses,  the  concentric 
cylinder  model  has  been  used  in  the  literature  due  to  its  simple  geometry.  Thermoelastic 
treatment  of  the  axisymmetric  concentric  cylinder  model  with  multiple  rings  are  given  in 
references  [16,  17].  An  elastic  analysis  of  a  multidirectional  coated  continuous  fiber 
composite  by  means  of  a  three-phase  concentric  cylinder  model  is  given  in  reference  [18]. 
Ebert  et  al.  [19]  and  Hecker  et  al.  [20]  were  the  first  to  introduce  elastic-plastic  behavior  for 
the  constituents  in  which  they  modeled  two  concentric  cylinders  and  verified  the  model  with 
experiments.  Gdoutos  et  al.  used  this  model  to  compute  thermal  expansion  coefficients 
[21  ]  and  stress-strain  curves  and  obtained  good  agreement  with  experimental  results  [22]. 
However,  this  approach  postulated  new  stress-strain  relations  in  which  the  strain  increments 
were  functions  of  the  stress  increment,  in  contrast  to  the  Pranatl-Reuss  relations  that  relate 
the  strain  increments  to  the  current  state  of  stress.  This  method  is  a  deformation  type  of 
theory  and  could  not  be  applied  to  nonproportional  loadings  and  ideally  plastic  material. 
More  recently  Lee  and  Allen  [23]  obtained  an  analytical  solution  for  elastic-perfectly  plastic 
fiber  and  matrix  obeying  Tresca's  yield  criterion. 

In  this  investigation  an  analytical  tool  was  developed  to  compute  the  three- 
dimensional  stress  state  in  a  composite  using  the  concentric  cylinder  model.  The  analysis 
accommodates  multiple  materials  having  elastic-plastic  behavior  with  strain  hardening.  The 
Prandtl-Reuss  (low  rule  with  the  von  Mises  yield  criterion  was  used.  In  addition, 
temperature-dependent  material  properties  were  taken  into  account.  The  analysis  was 
implemented  into  the  computer  code  FIDEP  (Finite  Difference  Code  for  Elastic-Plastic 
Analysis  of  Composites)  which  accounts  for  thermomechanica!  cyclic  loading. 


CHAPTER  3 
THEORY 


The  derivations  of  the  elastic-plastic  concentric  cylinder  equations  are  described  in 
this  section.  The  solution  technique  for  the  numerical  integration  of  the  Prandtl-Reuss 
equations  are  also  summarized. 

3.1  Concentric  Cylinder  Model 

A  representative  volume  element  of  the  composite  is  modeled  as  concentric  cylinders 
with  the  core  cylinder  representing  the  fiber  and  the  outer  ring  representing  the  matrix  (Fig. 
3.1).  The  fiber  and  matrix  radii  are  denoted  as  a  and  b,  respectively.  The  direction  of  the  z- 
axis  is  along  the  fibers  and  the  cylinders  are  infinitely  long  in  the  axial  direction. 
Cylindrical  coordinates  are  used  in  the  equations. 

The  following  assumptions  are  made  in  the  analysis.  The  temperature  distribution  is 
uniform  and  is  quasi-static.  A  perfect  bond  exists  between  the  constituents  of  the  composite 
so  that  there  is  no  slippage  or  separation  of  the  constituents.  The  concentric  cylinders  are  in 
generalized  plane  strain  and  are  subjected  to  axisymmetric  loadings  and  displacements  so 
that  the  shear  stresses  are  zero.  The  constituent  properties  are  isotropic.  The  fiber  is 
linearly  elastic.  The  matrix  follows  a  von  Mises  yield  surface  and  is  incompressible  in  the 
plastic  region,  i.e.,  hydrostatic  stresses  do  not  cause  plastic  deformation.  The  plastic 
deformation  of  the  matrix  is  governed  by  the  Prandtl-Reuss  flow  rule. 

The  following  boundary  conditions  are  imposed: 

1)  radial  stress  is  Pr  at  r=b, 

2)  finite  stresses  at  r={),  and 

3)  continuous  radial  displacements  and  radial  stresses  at  the  interface. 

3.2  Governing  Equations 

Equilibrium  and  compatibility  equations  in  cylindrical  coordinates  for  an 
axisymmetric  generalized-plane  strain  case  simplify  to  124]: 
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Fig.  3. 1  Concentric  Cylinder  Idealization  of  a  Unidirectional  Composite  and 
the  Concentric  Cylinder  Model. 


(1) 


dar  |  <7r  —  <7g  _  Q 
dr  r 

d£e  _  £r  ~  _  r> 

dr  r 


(2) 


and  the  stress-strain  equations  are  [25]: 

£r  =-=(°r-  +  <J,))  +  aT  +  £p  +  depr, 

E 

£g  =  ^(Gg-v(or  +  aJ)  +  aT  +  £p  +  d£p,  (3) 

E 

£,  =  ^f<72  -  V(<7r  +  Ge))  +  aT  +  £p  +  d£p 
k 


where  £rP,  EeP  and  £ZP  are  the  total  accumulated  plastic  strains  up  to,  but  not  including  the 
current  increment  of  loading,  d£rP,  d£@P  and  d£zP  are  the  plastic  strain  increments  due  to  the 
current  increment  of  loading,  £r,  £q  and  ez  are  the  total  strains,  ar,  cq  and  oz  are  the  stresses, 
a  is  the  secant  coefficient  of  thermal  expansion  (CTE),  v  is  the  Poisson’s  ratio,  E  is  the 
Young's  modulus  and  T  is  the  change  in  temperature  from  a  reference  state  in  which  the 
stresses  and  the  strains  are  assumed  to  be  zero. 

Substitution  of  Eqn.  3  into  Eqn.  2  to  eliminate  total  radial  and  hoop  strains  yields: 


+  W '<*r  +  Ge)-aET  +  E(£x  -  £p  -  depj)  +  ar  +  £p  +  d£p^j 


dr 

,  fl+vj  £p  +  d£pe  -  £p  -  d£p  _ 

■1  ~  ~  (  Gg  —  Gr)  +  - - e- - i - ^-  =  0 

Er  r 


(4) 


This  equation  together  with  the  equilibrium  equation  results  in  two  ordinary  differential 
equations  in  the  two  unknowns,  cr  and  Gq.  The  axial  strain,  Ez,  is  a  constant  value  across  the 
cross  section  because  of  the  imposed  generalized  plane  strain  condition. 
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3.3  Computation  of  the  Axial  Strain 

To  compute  the  axial  strain,  the  stress-strain  equation  in  the  axial  direction  (Eqn.  3) 
is  multiplied  by  E  r  and  integrated  over  the  cross  section: 

b  b  b  b  b 

£z  |  Erdr  -  j  a /dr  -  J  v(  oT  +  ae  )rdr  +  J  aETrdr  +J  E(  ef  +  dez  )rdr .  (5) 

0  0  0  0  0 

For  k  concentric  cylinders,  let  a;  be  the  outer  radius  of  the  ith  ring  and  aQ  =  0,  then  the 
integral  on  the  left  hand  side  reduces  to: 

j  Erdr  =  £2'2j~(  af  ~  4-i )  •  (6> 

0  i=l  ^ 

The  first  integral  on  the  right  hand  side  is  evaluated  using  the  global  equilibrium  equation  in 
the  axial  direction;  i.  e.  internal  forces  are  equal  to  the  external  applied  forces: 

b 

jo,(2nr)dr  =  Pt(nr2) ,  (7) 

o 

where  Pz  is  the  applied  axial  stress. 

Rearranging  the  terms,  the  equilibrium  equation  (Eqn.  1)  can  also  be  written  as  [26]: 

(ar  +  Ge)r  =  ^-(r2(Tr).  (8) 

dr 

Using  Eqn.  8,  the  second  integral  on  the  right  hand  side  of  Eqn.  5  becomes: 

b  d  * 

J  v ^ —  ( r 2 a r)dr  =  £  v,  f a2 or (a, )  -  a2. , o,  ( a, , )) ,  (9) 

o  ar  1=1 

Expanding  and  recollecting  terms  in  Eqn.  9: 
r  d 

J  v — ( r 2 a ,)dr  =  ]T  ( v, .  -  v(_ ,  )a2 a T ( at )  +  vkb 2 or ( b ) ,  (10) 

o  1=1 
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where  ar(b)  is  the  applied  radial  stress  at  r=b. 

The  remaining  terms  are  similarly  evaluated  assuming  constant  material  properties 
and  temperature  distribution  in  each  cylinder.  The  axial  strain  for  k  concentric  cylinders 
then  becomes: 


=  7- yi  ~2'LViCJr(al)(  K  ~  vM)-2vkor(b)+acEeT 
^ c  l  i= 1 

+Ty'LEi\<£*+de*)rdr  } 

0  i=\  0  ) 


(11) 


where,  Vj  is  the  volume  fraction  of  the  ith  cylinder,  Ec  is  the  axial  composite  modulus  and 
<Xc  is  the  axial  coefficient  of  thermal  expansion  of  the  composite  defined  by: 


and 


t 


In  the  case  of  two  concentric  cylinders  representing  an  elastic  fiber  and  an  elastic-plastic 
matrix,  k  =  2  and  ezP  =  0  in  the  fiber,  so  that  the  expression  for  the  axial  strain  simplifies  to: 


±-\pt~2V,<v,-  vm)o,(a)-2v„P,  +  a,E,T  +  —=■ 


j(£p  +dep)rdr  1,(12) 


where 


a'  =  Y(CCfEfVf  +  0CmEnVn)' 


Ec  =  EfVf  +  EmVm  , 
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and  Vf  and  Vm  are  the  fiber  and  matrix  volume  fractions,  respectively,  and  Pr  is  the  applied 
radial  stress  at  r=b.  In  Eqs.  11  and  12,  the  axial  strain  is  written  in  terms  of  the  radial 
stresses  at  the  interfaces.  The  total  plastic  strains  are  known  from  the  previous  loading 
steps  and  the  new  plastic  strain  increments  are  related  to  the  stresses  through  the  Prandtl- 
Reuss  flow  rule. 

Note  that  if  the  Poisson's  ratio  for  the  cylinders  is  the  same  and  the  applied  radial  stress  is 
zero,  the  second  and  third  terms  vanish  leaving  the  total  axial  strain  as  the  sum  of  the 
composite  elastic  strain  (e^),  composite  thermal  strain  (ezlhc)  and  composite  plastic  strain 
(e-/**),  i.  e.: 


£c  =  £ec  +  e'1*  +  £fe  =  —  +  acT  +  -4  —  f  (£p.  +  d£pt  )rdr . 

i  z  z  z  p  c  ,2  P  J  '  j  ‘  ' 


bl  £„ 


If  the  plastic  strain  is  constant  throughout  the  matrix,  then  we  have  for  the  total  strain: 


P  E  V 

=  —  +  a.T  +  -£-£■  ef"  , 

Ec  Ec 


where  £pm  is  the  average  matrix  plastic  strain. 

3.4  Finite  Difference  Formulation 

The  method  of  finite  differences  is  used  to  solve  the  two  ordinary  differential 
equations  [25].  The  disk  radius  is  divided  into  N  intervals  as  shown  in  Fig.  3.2.  There  are 
thus  N+l  stations,  the  first  station  being  at  the  center  of  the  disk  and  the  last  station  at  the 
outer  radius.  Eqs.  1  and  4  are  written  in  finite  difference  form  at  midpoints  of  these 
intervals  as  follows: 


r.  „  = 


o,.-oT 


r.  -  r. 


-  o 


r,  .-1 


2 


d  f  o$  ^  _  <V,/E, 

dr\E  ),_x  r-r  ., 
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r(k),  E(k),  a(k),  oy(k),  Ep(k) 
Ojj(k),  ejj(k),  AejjP(k) 


1 

1 

kA=INT 

1 

1 

■ 

kB=INT+1 

1 

1 

1 

Fiber/Mat  rix^^l 

1 

1 

1 

1 

1 

Interface  1 

1 

1 

1 

1 

Fig.  3.2 

Discretization  of  the  Concentric  Cylinder  Model. 

a+v)?JL  -  1  n+v^/fi-fl+v^.M,  ,_,/£„,  eic 
Er  i-'A  2  rt  -  ril/2 

In  this  manner  Eqs.  1  and  4  can  be  written  at  the  midpoints  of  n  stations  resulting  in  2n 
equations; 

/  =  2,...,/i  +  l,  (13) 

where 

=  (ri  ~  ri-x)l  2r,  ,  b,  =  E,  /  ,  c;  =  ri  lr , 

and 

A  =  fl  +  v,X v  +  a,)  , 

A=A+ v.Xl- v.  +a,), 

Ci  =  6in+vi-1Xvi.1-aicJ, 
tf,=6,(l+  v'_1)(l-  v,.,  -a.c,), 

A  =  £,T ( a,(l  +  v, )  -  a,„,  (1  +  v,_, )) , 

Pi  =  Ei(P-  +et(vi-vi_J), 

p;  =  +  <)-  ViJe!U  +  *£-.> 

+a,[<;  +  -  e£  (l/a,  + 1)  +  (1/a,  -  c, )] , 

where  the  plastic  strains,  £PJ ,  etc.,  are  the  updated  plastic  strains,  i.  e.  £PJ  =  ef,  +  dept,  etc. 

The  coefficients  at  the  left  hand  side  are  functions  of  material  properties  at  the  ith  station. 
Only  the  P*  term  on  the  right  hand  side  involves  plastic  strains.  The  unknowns  are  ar,j  and 
Oo,i,  i  =  l,...,n+l.  Using  the  boundary  conditions,  crn+i  =  Pr  and  oej  =  ar>i,  the 


—  1  "  °r,  +  +  a,°e,  -  0 

Ci 

GPri- 1  +  DPr,  +  ~  Fi°0.i  =  Qi-Pi 


11 


unknowns  reduce  to  2n;  orj,  i=l,...,n,  and  09,1,  i=2,...,n+l.  The  axial  stress  distribution,  oZ-l, 
i=l,...,n+l,  is  given  by  Eqn.  3. 


Moving  the  radial  stress  terms  in  Eqn.  1 1  to  the  left  hand  side  of  Eqn.  13  and  letting 

the  radial  stress  at  the  jth  interface  correspond  to  the  radial  stress  at  node  i=Ij,  i.  e., 
or(aj)  =  or  l  ,  the  second  expression  in  Eqs.  13  becomes: 


E  4-1 

+  D,Gr  i  +  H,a0._x  -  Fl(56  l  -  2(  v,  -  -  vl  +l) 


;= i 


,  (14) 


=  Q,  -  Ef  -  -  E,(  v,  -  vi_,)e! 


where 


1  k  r 

£,*  =  T(P,~2vkPr  +  cC'EJ  +  llb^E^e^rdr).  (15) 

Ec  J= i  o 

In  the  case  cf  two  concentric  cylinders,  k=2  and  j=l .  Let  INT  be  defined  as  the  index  of  the 
highest-numbered  node  in  the  fiber  and  INT+1  be  the  index  of  the  lowest-numbered  node 
in  the  matrix.  The  fiber-matrix  interface  lies  between  these  two  nodes.  Then,  for  all 
equations  in  which  i  *  INT+1,  i.  e.,  equations  for  nodes  not  adjacent  to  the  interface,  ni-nj.i 
is  zero  and  the  axial  strain  vanishes  from  these  equations.  For  the  INT+lst  equation,  Eqn. 
14  reduces  to: 


(2VfEm  I  Ec(vf  vm)  G/yr+i 

+  H m+\Ge,iNT  +  ^W+i ^e.wr+i  •  (16) 

=  Q/nt* i  ~  EmPfST+ 1  -*  Em(vm  —  Vf)£t 

One  final  step  before  Eqs.  13  are  written  in  matrix  form  is  to  eliminate  singular  terms  for 
i=2.  Note  that  Grji  =  Gej.  Thus  ori  and  09,1  terms  can  be  combined  in  Eqn.  1 3  to  obtain: 

(1  /  c2  +  a2)or ,  -  or  2  +  a2ce  2  =  0 
( H2  —  G2)ori  +  D2or2  —  F2gB2  —  Q2  —  P2 

In  the  first  equation  l/c'2  =  =  0.  In  the  second  equation,  the  term  H2-G2  is  expanded 

and  simplified.  The  final  equations  for  i=2  become: 
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a^OTy  Or  2  +  ,2  ~  ® 

a  +  *'•  )fl  -  2  v,  )h2cr, ,  +  D2CTr>2  -  F2ag2  =  Q2-P2 


Hence  Eqn.  13  is  written  in  matrix  form  as: 

A  x  =  B  ,  (17) 

where  A  is  a  2n  by  2n  matrix  of  constant  coefficients,  x  is  the  radial  and  hoop  stress  vector 
of  length  2n  and  B  is  a  vector  of  length  2n,  as  shown  in  Fig.  3.3.  In  matrices  A  and  B  all 
the  constants  are  known  except  the  plastic  strain  increments  which  are  presumed  or 
computed  from  the  Prandtl-Reuss  relations. 

3.5  Plasticity  Theory 

The  plastic  strains  increments  are  computed  using  Prandtl-Reuss  relations  [25].  In 
this  section  two  forms  of  the  Prandtl-Reuss  equations  are  summarized;  equations  relating 
plastic  strain  increments  to  stresses  and  equations  relating  plastic  strain  increments  to 
modified  total  strains. 

3.5.1  Prandtl-Reuss  Relations 

Prandtl  and  Reuss  assumed  that  the  plastic  strain  increment  is  proportional  to  the 
instantaneous  stress  deviation,  i.  e.: 

=  (18) 

where  Sij  is  the  deviatoric  stress  tensor  and  X  is  a  nonnegative  constant  which  may  vary 
throughout  the  loading  history.  These  equations  imply  that  the  plastic  strain  increments 
depend  on  the  current  stress  state  and  not  on  the  stress  increment  required  to  reach  this 
state.  To  determine  X,  both  sides  are  squared  and  multiplied  by  2/3  to  obtain: 

o’) 

Define  effective  plastic  strain  increment  and  effective  stress,  respectively,  as: 
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Figure  3.3  The  Linear  System  of  Equations  Resulting  from  the  Finite  Differences  Formulation  of  the 
Governing  Equations  for  the  Concentric  Cylinder  Model. 


*'  =  ^ §<*£  and  <r«  =  ffiSe.  (20) 

Using  these  definitions  X  is  determined  from  Eqn.  20  as: 

A=-—  (21) 

2o<ff 

and  the  Prandtl-Reuss  relations  become: 

3 

=  — -—S  .  (22) 

These  equations  are  used  together  with  the  von  Mises  yield  criterion.  Yielding  begins  when 
the  effective  stress  reaches  the  yield  stress  determined  from  a  uniaxial  tensile  test. 

3.5.2  Plastic  Strain-Total  Strain  Plasticity  Relations 

The  above  form  of  the  Prandtl-Reuss  relations  relate  the  plastic  strain  increments  to 
the  stresses.  Another  method  to  computing  the  plastic  strain  increments  is  to  use  equations 
that  relate  the  plastic  strain  increments  to  the  modified  total  strains.  These  equations  will  be 
more  stable  with  respect  to  the  loading  increment  size  and  converge  faster  for  most  loading 
cases.  In  this  approach  the  plastic  strain  increments  are  computed  using  modified  Prandtl- 
Reuss  relations  derived  by  Mendelson  [25]. 

Assume  a  loading  path  to  a  given  state  of  stress  and  total  plastic  strains  EijP  and  let 
the  next  load  step  be  applied  producing  additional  plastic  strains  AEjjP.  The  total  strains  can 
be  written  as: 


=  e*  +  <£  +  <  +  te-r  (23) 

where  £jjc  is  the  elastic  component  of  the  total  strain,  Ejjlh  is  the  thermal  strain,  £jjP  is  the 
accumulated  plastic  strain  up  to  (but  not  including)  the  current  increment  of  load,  and  AEjjP 
is  the  increment  of  plastic  strain  due  to  the  increment  of  load.  The  previous  plastic  strains 
EjjP  are  presumed  to  be  known,  and  AEjjP  is  to  be  calculated. 


15 


The  modified  total  strains  are  defined  as: 


(24) 


The  deviatoric  strains  are  obtained  by  subtracting  the  mean  strain  trom  the  diagonal 
compo'  ent  of  both  sides: 

=  <  +  (25) 


where  £jje  is  the  elastic  strain  deviator  tensor  and  £jj’  is  the  modified  strain  deviator  tensor. 
From  Hooke's  Law  and  Prandtl-Reuss  Eqn.  1 8: 


Ml 

2 GX  ' 


(26) 


where  G  is  the  shear  modulus.  Using  this  expression  and  Eqn.  18,  Eqn.  25  becomes: 


e'  =  (1  +  — —  )dep  . 
"  V  2GA  “ 


Multiplying  both  sides  of  the  equation  by  two  thirds  itself: 


-e'  „  e’  „  =  -(1  + 

3  '  3  2GA 


*/  V 


Define  effective  modified  total  strain  in  a  similar  fashion  to  the  effective  incremental  plastic 
strain  rate  to  obtain: 


•s 


=  (1  + 


2GA 


)dep. 


The  modified  Prandtl-Reuss  relations  then  become: 


de 


p 

V 


'27) 


(28) 
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where  eij'  are  the  modified  total  deviatoric  strains,  e'cff  is  the  equivalent  or  effective  modified 
total  strain  defined  by: 


iff- 


(29) 


and  dep  is  the  equivalent  or  effective  plastic  strain  increment  defined  by  Eqn.  20.  The 
modified  Prandtl-Reuss  equations  in  expanded  form  are: 


dep 

d£px=TT-(2£\-e'y-e\), 

dep 

drp  =  — —(2e'  -£'  -£'  ), 

y  1  y  1  * h 

■>fc  iff 

dep 

d£p  =  -f-(2e\  -£\  -£'y )  =  -(de[  +  d£p) 

i£  iff 


d£p  d£p 

dep  =^£j£-£'  t  dep  =  “eJLe’  t 

iff  iff 


dep  =f£ ‘Jl£- 
u  3e' 

iff 


Equation  28  is  equivalent  to  the  Prandtl-Reuss  Eqn.  18,  but  relate  the  incremental  plastic 
strains  to  modified  total  strains  instead  of  the  stresses. 


The  relationship  between  the  effective  incremental  plastic  strain  and  effective 
modified  total  strain  is  obtained  by  substituting  \  from  Eqn.  21  in  Eqn.  27: 


e  iff  -  d£<ff  +  a't  ■ 


(30) 


Let  the  effective  stress  from  the  previous  loading  step  be  oeJ) .  By  expanding  the  effective 
stress  in  a  Taylor  series  about  atff : 


aui  = 


do 


iff 


d£p 

utiff 


d£p 


A~1 


the  effective  stress  is  eliminated  from  Eqn.  30  to  obtain: 
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(31) 


d£'ff  ~ 


£‘ff  3^ 


The  effective  modified  total  strain  has  no  physical  significance  and  is  purely  a  mathematical 
quantity  even  in  the  uniaxial  case.  In  Eqn.  31  a  small  error  in  effective  modified  total  strain 
or  the  effective  stress  will  produce  the  same  order-of-magnitude  error  in  the  effective  plastic 
strain  increment  since  the  denominator  is  approximately  unity.  Equations  28,  29  and  31  are 
used  simultaneously  to  determine  plastic  strain  increments  at  each  loading  step. 

3.6  Solution  Technique 

Solution  of  the  elastic-plastic  concentric  cylinder  problem  requires  the  solution  of  5 
equations;  the  equilibrium  and  compatibility  Eqs.  1  and  6,  and  the  three  Prandtl-Reuss 
relations  (Eqn.  22  or  28).  After  elimination  of  the  strains  using  the  constitutive  relations 
there  are  6  unknowns  left;  or,  <Jq,  oz,  derp,  deep,  d£zp.  The  axial  stress  is  eliminated  using 
Eqn.  1 1  reducing  the  number  of  unknowns  to  5.  The  ordinary  differential  equations  (Eqs.  1 
and  6)  are  solved  using  the  method  of  central  finite  differences  to  obtain  a  linear  system  of 
equations  (Eqn.  13).  These  simultaneous  equations  are  used  in  an  iterative  plasticity 
algorithm  to  solve  for  the  stresses. 

There  are  two  iterative  algorithms  that  were  considered  for  this  investigation  both  of 
which  use  forward-Euler  integration  scheme.  The  first  algorithm  uses  the  Prandtl-Reuss 
relations  (Fig.  3.4).  For  an  increment  of  load,  a  distribution  is  assumed  for  the  plastic  strain 
increments.  The  finite  difference  equations  are  then  solved  for  a  first  elastic  approximation 
to  the  stresses  and  strains.  At  the  same  time,  the  effective  plastic  strain  increment  is 
calculated  (Eqn.  20).  From  the  stress-strain  curve  the  corresponding  effective  stress  is 
computed.  Finally,  a  new  estimate  computed  for  the  p las  -  -t-jm  increments  and  these 
steps  are  repeated  until  convergence  is  obtained.  T'r.  ■  •  m  is  not  a  very  stable 

algorithm  and  may  not  converge  for  large  load  increments. 

The  second  iterative  method  uses  the  modified  version  of  the  Prandtl-Reuss  Eqn.  28 
or  the  plastic  strain-total  strain  relations  developed  by  Mendelson  (251.  This  method  is  also 
called  the  elastic-predictor  radial  corrector  method.  In  this  algorithm  (Fig.  3.5),  total  strains 
are  obtained  from  the  stresses  from  Eqn.  13  with  assumed  or  previously  calculated  plastic 
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Figure  3.4  Schematic  of  the  Iterative  Technique  for  Computing  Plastic  Strains  Using  the 
Prandtl-Reuss  Relations. 
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Figure  3.5  Schematic  of  the  Iterative  Technique  for  Computing  Plastic  Strains  Using  the 
Modified  Prandtl-Reuss  Relations. 
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strain  increment  estimates.  Modified  total  strains  and  effective  total  strain  are  obtained  from 
Eqs.  38  and  stress-strain  curve  using  Eqn.  43.  The  new  plastic  strain  increments  are  then 
determined  from  Eqn.  40.  This  process  is  repeated  until  convergence  of  the  plastic  strains 
is  obtained.  This  method  converges  for  even  large  loading  increments  and  converges  very 
rapidly. 


Both  of  these  algorithms  were  implemented  into  the  program  FIDEP,  Finite- 
Difference  Code  for  Elastic-Plastic  Analysis  of  Composites.  However,  since  the  first 
algorithm  did  not  offer  any  advantage  and  was  unstable  except  for  very  small  time  steps,  the 
method  using  the  modified  plastic-strain  total-strain  relations  was  used. 
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CHAPTER  4 

IMPLEMENTATION  IN  FIDEP 


4.1  Algorithm 

The  algorithm  and  the  pseudo-code  for  computing  the  stresses  and  plastic  strain 
increments  using  the  modified  strain  theory  as  implemented  into  the  FORTRAN  code 
FIDEP  are  shown  in  Figures  4.1  and  4.2.  The  procedure  for  calculating  the  stresses  and 
strains  is  as  follows: 

1 .  A  stress  increment  and  a  temperature  increment  are  applied. 

2.  The  finite  difference  equations  (13)  are  solved  for  the  new  stress  state 
assuming  an  elastic  deformation,  i.  e.,  derP  =  deeP  =  dezP  =  0. 

3.  A  new  yield  surface  is  computed  for  the  present  temperature  at  all  stations. 

4.  Effective  stress  is  compared  to  this  new  yield  surface  for  each  station. 

5.  If  the  effective  stress  is  less  than  the  effective  stress  at  all  nodes  the  plasticity 
subroutine  is  bypassed  to  go  to  the  next  thermal  and  mechanical  load 
increments.  Otherwise  a  nonzero  plastic  strain  increment  is  assumed. 

6.  The  linear  system  of  equations  (Eqn.  17),  with  the  assumed  plastic  strain 
increment,  are  solved  for  the  radial  and  hoop  stresses.  The  axial  stress  is 
computed  from  these  stresses. 

7.  The  modified  total  strains  and  the  equivalent  total  strain  are  computed  from 
Eqs.  24  and  29,  respectively. 

8.  The  effective  incremental  plastic  strain  is  computed  from  Eqn.  20  and  the 
new  incremental  plastic  strains  are  computed  from  the  plastic  strain-total 
strain  relations  (Eqn.  28). 

9.  These  strains  are  then  compared  with  the  previous  plastic  strain  increments 
and  if  the  difference  is  less  than  a  specified  tolerance,  go  to  step  1. 
Otherwise  the  linear  system  of  equations  (Eqs.  17)  are  solved  with  the  new 
incremental  plastic  strains. 

10.  Steps  6-9  are  repeated  until  convergence  of  the  incremental  plastic  strains  is 
obtained. 
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•  temperature  Dependent  Material  Properties 

•  Loading  History 
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Figure  4.1  Elastic-Plastic  Algorithm  Used  in  FIDEP. 
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Figure  4.2  Pseudo-Code  for  the  Computer  Program  FIDEP. 
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Figure  4.2  Pseudo-Code  for  the  Computer  Program  FIDEP  (Continued). 


4.2  Faster  Convergence  Scheme 


To  reduce  the  number  of  iterations,  a  faster  convergence  scheme  was  adopted  [24]. 
In  this  method  the  three  previous  plastic  strain  increment  estimates  and  their  differences  are 
saved  and  extrapolated  to  a  plastic  strain  increment  that  would  result  in  a  difference  of  zero. 
In  Fig.  4.3,  the  slope  of  the  line  is  calculated  and  the  new  strain  is  obtained  as  -b/a,  which  is 
the  y-i ntercept  of  the  plastic  strain  vs.  the  difference  graph. 


S  =  aep  +  b 

=  £*2;- 
S2=e’°>  -e  *2> 
b 

£  =  — 

a 

where 

5,  ~  ,  A 

e'  e*3*  -e*2*  ~ 1  d2 
b  =  S2 -a£p<v 


Figure  4.3  Schematic  for  the  Computation  of  a  New  Estimate  for  the  Plastic  Strains 
Using  the  Previous  Estimates  and  their  Differences. 
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CHAPTER  5 
PROGRAM 


5.1  General  Description 

FIDEP  (Finite-Difference  Code  for  Elastic-Plastic  Analysis)  is  a  user-friendly 
computer  program  designed  to  solve  the  problem  of  concentric  cylinders  subjected  to 
thermal  and  mechanical  loading.  The  program  is  simple  in  structure  and  scope  since  it  was 
intended  to  be  used  as  a  practical  research  tool  for  the  determination  of  stress  concentrations 
in  the  matrix  around  the  fiber  in  unidirectional  continuous  fiber  reinforced  metal  matrix 
composites.  The  features  of  the  program  include: 

•  generalized  plane  strain  and  axisymmetric  geometry, 

•  small  displacements, 

•  isotropic  and  incompressible  fiber  and  matrix, 

•  elastic  fiber  and  time-independent  bilinear  elastic-plastic  matrix 

(including  perfectly  plastic  matrix), 

•  perfect  bonding  between  the  fiber  and  the  matrix, 

•  uniform  temperature  change  throughout  the  composite. 

More  advanced  features  such  as  time-dependent  material  behavior,  multiple 
concentric  cylinders  are  absent  from  the  current  version  of  the  program  since  only  speedy 
prediction  of  the  time-independent  nonlinear  behavior  of  a  two-phase  unidirectional 
composite  was  required  at  the  time. 

Operation  of  FIDEP  program  is  straightforward,  requiring  the  modification  of  a 
loading  history  file  and  a  material  properties  data  file.  All  input  data  are  in  free  format  and 
some  descriptive  titles  are  allowed.  The  following  sections  describe  the  operation  of  the 
program. 

5.2  Program  Operation 


The  execution  of  the  program  consists  of  changing  two  input  files  and  executing  the  code. 
The  input  files  are  the  loading  history  file,  FDLOAD.DAT,  and  temperature-dependent 
constituent  property  data  file,  FDMAT.DAT.  Running  the  code  creates  four  output  files: 


FDOUTS,  FDOUTE,  FDUSP,  FDSUM.  The  first  two  files  consist  of  the  stress  and  strain 
history,  respectively,  of  the  matrix  at  the  interface.  The  output  file  FDUSP  consists  of  the 
stress-strain  distribution  that  exists  in  the  composite  at  the  end  of  the  imposed  loading 
history.  The  output  file  FDSUM  consists  of  a  summary  of  the  fiber  and  matrix  stresses  at 
the  interface  at  the  end  of  each  half-cycle.  In  the  following  sections  the  input  and  output  file 
formats  are  explained  in  detail. 

5.2.1  Input 

The  input  consists  of  two  files:  loading  file  and  material  properties  file.  The  loading 
file  data  specifies  the  number  of  steps,  temperature,  and  applied  mechanical  stresses  and  the 
second  file  specifies  the  temperature  dependent  mechanical  properties  for  the  constituents. 

The  loading  file,  FDLOAD,  format  is  listed  in  Table  5.1.  An  example  loading 
history  data  file  for  cyclic  loading  at  room  temperature  after  cooldown  from  processing  is 
shown  in  Fig.  5.1. 

The  material  property  data  specifies  the  following  properties  for  the  matrix  and  the 
fiber:  elastic  modulus,  coefficient  of  thermal  expansion,  Poisson's  ratio,  yield  stress  and 
plastic  modulus  for  each  temperature.  The  material  data  file,  FDMAT,  format  is  listed  in 
Table  5.2.  An  example  material  data  file  for  SCS-6  fiber  and  Ti-24A1-1  INb  matrix  material 
is  shown  in  Fig.  5.2. 

5.2.2  Output 

There  are  four  output  files:  FDOUTS.DAT,  FDOUTE.DAT,  FDUSP.DAT  and 
FDSUM.DAT.  The  first  two  files  contain  stresses  and  strains,  respectively,  in  the  matrix  at 
the  fiber/matrix  interface  at  each  step  increment.  FDUSP.DAT  contains  the  stresses  and 
strains  along  the  radius  for  the  final  loading  step.  FDSUM.DAT  contains  a  summary  of 
the  stresses  and  strains  at  the  interface  in  the  fiber  and  the  matrix  at  the  end  points  of  each 
cycle.  The  runtime  is  printed  on  the  screen  at  the  completion  of  the  program.  The  output 
formats  are  described  in  the  following  paragraphs. 

FDOUTS.DAT  echos  the  inputs  from  the  loading  and  material  data  files  at  the  beginning  of 
the  file.  Following  the  input,  the  stresses  and  the  mechanical  strain  in  the  matrix  computed 
at  the  fiber/matrix  interface  are  printed  in  the  following  order:  computational  step,  and 
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Cooldown  and  Mechanical  Cycling  at  Room  Temperature 
Comment  line  2 
Comment  line  3 
Comment  line  4 


o 


0 

30 

60 

90 

120 

0 

0 

0 

0 

1010  25 

25 

25 

25 

0 

0 

0 

0 

0 

0 

350 

35 

350 

0 

0 

0 

0 

Figure  5.1  Example  Loading  History  Input  File,  FDLOAD.DAT. 


TI-24-11  Matrix  and  SCS-6  Fiber  Properties 
Comment  line  2 
2 

11 


TEMP  E(GPa)  CTE(lE-6/C) 

MU  SY(MPa)  EP(GPa) 

20  413 

4.72 

0.3 

1.E6 

0. 

101  413 

4.81 

0.3 

1.E6 

0. 

203  413 

4.96 

0.3 

1.E6 

0. 

299  413 

5.09 

0.3 

1.E6 

0. 

400  413 

5.23 

0.3 

1.E6 

0. 

500  413 

5.36 

0.3 

1.E6 

0. 

598  413 

5.44 

0.3 

1.E6 

0. 

702  413 

5.51 

0.3 

1.E6 

0. 

800  413 

5.66 

0.3 

1.E6 

0. 

900  413 

5.73 

0.3 

1.E6 

0. 

1010  413 

Q 

5.75 

0.3 

1.E6 

0. 

TEMP  E(GPa)  CTE(lE-6/C) 

MU 

SY(MPa)  EP(GPa) 

20  84.10 

11.33 

0.3  950.0 

0.00 

315.  88.40 

11.88 

0.3  410.0 

3.29 

482.2 

68.26 

12.22  0.3 

333.3 

3.71 

649.0 

48.09 

12.73  0.3 

256.5 

4.12 

704.4 

42.08 

12.95  0.3 

214.1 

3.89 

760.0 

36.07 

13.22  0.3 

171.8 

3.67 

885.0 

23.66 

13.93  0.3 

112.8 

3.91 

1010.0 

11.25 

14.97  0.3 

53.8 

4.15 

Reference  Temperature 
1010 

Vf  b 
0.35  1.0 


Figure  5.2  Example  Material  Properties  Input  File,  FDMAT.DAT. 
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Table  5.1  Format  for  the  Loading  Data  File,  FDLOAD.DAT. 


Card 

Column 

Variable 

Definition 

14 

Comment  lines 

5 

1 

NLOAD 

Number  of  columns  of  loading  data  to  be  read 

6 

1-NLOAD 

NOT 

Step  number 

7 

1-NLOAD 

7E 

Temperatures  (°C)  corresponding  to  step  NOT 

8 

1-NLOAD 

PAT 

Applied  axial  stresses  (MPa)  corresponding  to  step  NOT 

Table  5.2  Format  for  the  Material  Properties  File,  FDMAT.DAT. 


Card 

Column 

Variable 

Definition 

1-2 

Comment  lines 

3 

1 

NOMAT 

Number  of  constituents  ("* ) 

4 

1 

NROW 

Number  of  rows  of  material  property  data  for  the  first 

constituent  following  header  line 

5 

1 

Header  line  for  the  properties 

6 

1 

TEMP 

Temperature  (°C) 

2 

EMAT 

Modulus  (GPa) 

3 

CTEMAT 

Coefficient  of  thermal  expansion  (E-6/°C)  (2) 

4 

MUMAT 

Poisson's  ratio 

5 

SYMAT 

Yield  stress  (MPa)  (3) 

6 

EPLMAT 

Plastic  modulus  (GPa)  (4) 

7- 

Repeat  cards  4-6  for  the  matrix 

8 

Comment  line 

9 

1 

TREF 

Reference  temperature  (5) 

10 

Comment  line 

11 

1 

VF 

Fiber  volume  fraction 

2 

B 

Outer  radius  of  matrix  (usually  1 .0) 

Notes:  (1)  The  number  of  constituents  is  limited  to  two  at  the  present. 

(2)  Secant  coefficient  of  thermal  expansion  (CTE)  is  used.  The  reference  temperature  for 
the  secant  CTE  should  be  the  initial  processing  temperature. 

(3)  (4)  For  the  fiber,  artificially  high  value  is  given  tor  the  yield  stress  and  zero  value  is 
given  for  the  plastic  modulus  since  the  fiber  is  assumed  to  be  elastic. 

(5)  This  is  the  same  temperature  as  the  initial  processing  temperature  used  in  loading 
history  data  file. 
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corresponding  temperature,  applied  axial  stress,  effective  stress,  radial  stress,  hoop  stress, 
axial  stress,  original  yield  surface,  and  axial  strain.  The  column  titles  are  shown  below. 

Step  Temp  Oapp  Oe(f  0r  O0  Oz  Y.  S.  £z 

The  stresses  are  given  in  units  of  MPa,  temperature  in  °C,  and  the  strain  in  mm/mm.  The 
effective  stress  is  defined  by: 

aeff  =  /  2[(ar  -  cr^ +  (or  -  a,  )2  +  {ot  -  ogf] 

which  is  equivalent  to  Eqn.  20. 

FDOUTE.DAT  consists  of  the  strains  in  the  matrix  at  the  fiber/matrix  interface  in  the 
following  order:  computational  step,  and  corresponding  temperature,  applied  axial  load, 
radial  plastic  strain,  hoop  plastic  strain,  axial  plastic  strain,  radial  mechanical  strain,  hoop 
mechanical  strain,  axial  mechanical  strain  and  thermal  strain.  The  column  titles  are  shown 
below: 


Step  Temp  Sapp  £r  £z  ^rme  ^Gme  Ezme  eth 

FDUSP.DAT  consists  of  the  stresses  and  plastic  strains  across  the  cross  section  of  the 
composite  at  the  end  of  loading  history  in  the  following  order:  radius,  effective  stress,  radial 
stress,  hoop  stress,  axial  stress,  radial  plastic  strain,  hoop  plastic  strain  and  axial  plastic 
strain.  The  column  titles  are  shown  below: 

Radius  CJcff  Or  Oq  Oz  Ejp  £©p  E^p 

EDSUM.DAT  consists  of  a  summary  of  the  stresses  and  strains  at  the  interface  in  the  fiber 
and  matrix  at  the  end  points  of  each  loading/unloading  cycle  in  the  following  order:  step 
number  corresponding  to  steps  in  the  loading  file,  FDLOAD,  temperature,  applied  stress, 
axial  fiber  stress,  axial  matrix  stress,  effective  matrix  stress,  effective  matrix  strain  and  axial 
mechanical  strain  in  the  matrix.  The  column  titles  are  shown  below: 

Step  Temp  ®app  ®eff  ^zmo 
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5.2.3  Program 


The  program  FIDEP  is  listed  in  Appendix  A.  The  program  uses  a  default  value  of  50 
divisions  in  the  cylinders  for  applying  the  finite  difference  method.  These  finite  difference 
equations  result  in  a  linear  system  of  100  equations  in  100  unknowns  represented  by  the 
variable  NRA  in  the  program.  If  a  better  accuracy  is  desired  NRA  can  be  changed  to  200  in 
the  subroutine  SOLVE. 

5.2.4  Execution  of  FIDEP  on  a  VAX/VMS  Machine 

The  program  FIDEP  can  be  run  interactively  or  in  a  batch  mode  on  a  VAX-VMS  machine 
and  requires  IMSL  library  of  subroutines  [27].  Before  execution,  the  program  is  compiled 
and  linked  on  a  VAX/VMS  machine  as  follows: 

>  FOR  FIDEP 

>  LINK  FIDEP+IMSL/LIB 
To  run  FIDEP  in  real  time  type 

>  RUN  FIDEP 

To  run  FIDEP  on  a  batch  mode  a  command  file  has  to  be  created.  In  this  command  file  the 
input  and  output  files  can  be  renamed.  An  example  command  file  is  shown  in  Fig.  5.3. 
This  command  file  is  then  submitted  to  the  batch  using  the  command 

>  SUBMIT/NOPRINT  filename 

After  the  job  has  been  completed  a  log  file  will  appear  in  the  main  directory. 


^  ^  ^  ^ 


SET  VERIFY 

ASSIGN/USER.MODE  [COKERDU.FIDEP] FDMAT.DAT  FDMAT 
ASSIGN/USER_MODE  [COKERDU.FIDEP3FDLOAD.DAT  FDLOAD 
ASSIGN/USER.MODE  [COKERDU.FIDEP3FDOUTS.DAT  FOOUTS 
ASSIGN/USER.MODE  [COKERDU.FIDEP] FDOUTE.DAT  FDOUTE 
ASSIGN/USER_MODE  [COKERDU . FIDEP]FDUSP. DAT  FDUSP 
ASSIGN/USER.MODE  [COKERDU.FIDEP3FDSUM.DAT  FDSUM 
RUN  [COKERDU . FIDEP]FIDEP3B 


Figure  5.3  Example  Command  File  for  Running  Batch  Jobs  on  a  VAX/VMS  Machine. 


CHAPTER  6 

DEMONSTRATION  PROBLEMS 


Several  problems  are  solved  using  FIDEP  code  to  illustrate  the  capability  for 
computing  micromechanical  stresses  and  to  compare  with  solutions  obtained  by  other 
methods.  The  first  case  is  the  comparison  of  FIDEP  solution  to  the  closed  form  solution 
for  an  elastic  problem.  The  second  case  is  the  comparison  of  FIDEP  results  with  finite 
dement  analysis  results  for  thermal  cool-down  and  thermal  cyclic  loading.  The  third  case  is 
an  application  problem  in  which  the  stresses  and  strains  are  predicted  in  a  unidirectional 
composite  subjected  to  in-phase  and  out-of-phase  thermomechanical  fatigue  behavior  and 
the  strains  are  compared  with  experimental  measurements. 

6.1  Material  Properties 

In  these  computations  two-material  concentric  cylinder  model  consisting  of  SCS-6 
silicon  carbide  fiber  and  Ti-24A1-1  INb  matrix  is  used.  The  mechanical  properties  for  SCS- 
6  and  Ti-24A1- 1  INb  are  shown  in  Table  6.1  as  a  function  of  temperature.  The  SCS-6  fiber 
is  isotropic  and  elastic  with  only  the  coefficient  of  thermal  expansion  varying  with 
temperature.  An  adequate  representation  of  the  Ti-24A1-1  INb  matrix  was  attained  using  a 
bilinear  elastic-plastic  model  with  temperature  dependent  elastic  and  plastic  moduli,  yield 
stress  and  coefficient  of  thermal  expansion. 

6.2  Comparison  with  an  Elastic  Solution 

An  elastic  closed  form  solution  for  two  concentric  cylinders  and  generalized 
axisymmetric  condition  is  presented  in  Appendix  2  [28].  The  stresses  at  the  cross-section 
obtained  by  FIDEP  and  by  the  closed  form  solution  for  AT=-100°C  is  shown  in  Fig.  6.1. 
In  this  analysis  temperature  independent  elastic  properties  were  used  with  the  following 
properties:  Ef  =  414  GPa,  Vf  =  0.22,  ctf  =  4.7E-6  f° C,  Em  =  90  GPa,  vm  =  0.30,  am  =  10.7E- 
6  /°C,  where  f  denotes  the  fiber  and  m  denotes  the  matrix.  The  fiber  volume  fraction  was 
35%.  The  stresses  obtained  by  the  two  method  are  in  excellent  agreement  as  shown  in  Fig. 
6.1. 


Stresses  at  the  cross-section  using  the  elastic-plastic  and  the  elastic  solution  is 
compared  in  Fig.  6.2  using  the  mechanical  properties  from  Table  6.1.  The  temperature 
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Table  6.1 


Mechanical  Properties  for  SCS-6  Fiber  and  Ti-24AI-llNb  Matrix. 


SCS-6  Fiber  Properties 
v  =  0.22 
E  =  414  GPa 


Temperature  a  *  (10'6/°C) 
(°C) 


20 

4.70 

93 

4.81 

204 

4.97 

316 

5.12 

427 

5.26 

538 

5.38 

649 

5.50 

760 

5.60 

871 

5.70 

982 

5.78 

1010 

5.80 

Ti-24AI-1 1  Nb  Matrix  Properties 
v  =  0.22 


Temperature 

(°C) 

a  *  (10'6/°C) 

Elastic 

Modulus 

(GPa) 

Yield  Stress 
(MPa) 

20 

12.33 

94 

604 

93 

12.47 

92 

560 

204 

12.78 

91 

498 

316 

13.21 

89 

447 

427 

13.75 

79 

421 

538 

14.42 

70 

381 

649 

15.20 

50 

357 

760 

16.10 

25 

252 

871 

17.11 

18 

138 

982 

18.25 

16 

38 

1010 

18.55 

15 

30 

Plastic 

Modulus 

(GPa) 


1.3 
0.9 
0.7 
0.7 
0.4 
0.1 

0 

2.3 
2.6 
1.2 
1.0 


*  a  is  secant  coefficient  of  thermal  expansion  with  a  reference  temperature  of  1010°C. 


Stress  (MPa) 


Radial  Distance  from  the  Fiber  Center 


Figure  6.1  Variation  of  Elastic  Stresses  with  Radius  Obtained  by  Elastic  Closed-Form  Solution 
and  FIDEP  for  AT  =  -  100°C. 
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Stress  (MPa) 


Nondimensional  Radius  (r/b) 

Figure  6.2  Variation  of  the  Stresses  Across  the  Cross-Section  of  the  CCM  After  Cool-Down 
from  Processing  Temperature  of  1010°C  for  Elastic  and  Elastic-Plastic  Matrix. 
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dependent  elastic  solution  was  obtained  from  FIDEP  by  making  the  yield  stress  artificially 
high  (1.E6  MPa)  for  the  matrix.  The  loading  input  file  and  a  schematic  of  the  loading 
history  is  shown  in  Fig.  6.3.  The  yield  surface  is  first  reached  at  the  interface  and  the 
matrix  near  the  interface  becomes  plastic.  As  the  thermal  loading  increase  the  plastic  zone 
spreads  outwards  to  the  rest  of  the  matrix  and  a  redistribution  of  the  stresses  occur.  In  the 
elastic  region,  the  radial  and  hoop  c messes  change  proportional  to  1/r2  and  the  axial  stress 
remain  uniform.  Inside  the  plastic  zone,  the  radial  and  hoop  stresses  remain  almost  uniform 
and  the  axial  stress  decreases  when  the  effective  stress  cannot  extend  beyond  the  yield 
stress.  The  remaining  load  is  redistributed  to  the  fiber  and  the  elastic  region  of  the  matrix. 

6.3  Comparison  with  Finite  Element  Method 

Two  examples  are  presented  and  the  results  from  FIDEP  are  compared  with  finite 
element  method  (FEM)  results  obtained  from  nonlinear  finite  element  code  MAGNA  [29]. 
The  temperature  dependent  plasticity  subroutines  were  implemented  in  MAGNA  by  J.  L. 
Kroupa  [30].  The  geometry  and  the  boundary  conditions  for  the  finite  element  model  of  an 
axisymmetric  concentric  cylinder  geometry  under  generalized  plane  strain  condition  are 
shown  in  Fig.  6.4  [30].  Contact  elements  were  used  to  model  the  fiber-matrix  interface, 
which  allowed  the  separation  of  the  fiber  and  matrix  when  the  matrix  stresses  at  the  interface 
became  positive.  Isotropic  hardening  model  with  a  bilinear  elastic-plastic  stress-strain  curve 
was  used  to  model  the  matrix.  The  fiber  was  elastic  and  the  volume  fraction  was  35%.  The 
material  properties  listed  in  Table  6. 1  were  used  for  the  analysis. 

6.3.1  Cool-Down  from  Processing  Temperature 

The  first  comparison  problem  with  FEM  was  the  thermal  cool-down  associated  with 
processing  the  material.  The  loading  input  file  and  a  schematic  of  the  loading  history  is 
shown  in  Figure  6.3.  The  fiber  and  matrix  are  assumed  to  have  zero  stresses  and  strains  at 
the  initial  consolidation  temperature  of  1010°C.  As  the  material  cools  down,  the  modulus 
changes  with  temperature  and  the  difference  in  coefficient  of  thermal  expansion  between  the 
fiber  and  matrix  leads  to  thermal  stresses  in  the  fiber  and  matrix.  The  material  properties 
are  those  shown  in  Table  6.1.  The  output  from  the  program  is  presented  in  Appendix  C. 
The  files  FDOUTS.DAT  and  FDUSP.DAT  were  used  in  Figs.  6.5-6.8  discussed  below. 

The  three-dimensional  stresses  computed  during  the  cool-down  process  are 
illustrated  in  Fig.  6.5.  The  computed  stresses,  shown  as  the  dashed  lines  with  symbols  in 


Figure  6.3 


Loading  File  1 

Cool-Down  to  Room  Temperature 

From  a  Processing  Temperature  of  1010°C 

For  SCS-6/Ti-24Al-llNb  Unidirectional  Composite 

2 

0  30 

1010  20 

0  0 


(a) 


Temperature  (°C) 


(b) 


a)  Loading  Input  File,  and  b)  Schematic  of  the  Temperature  History  for  a 
Thermal  Cool-Down  Simulation  of  SCS-6/Ti-24Al-l  INb  Unidirectional 
Composite. 
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v  -  uniform 


V7777777777Z^77777, 

v  =  0 


Fig.  6.4  Finite  Element  Model  and  the  Boundary  Conditions  for  an  Axisymmetric 
Concentric  Cylinder  Geometry  Under  Generalized  Plane  Strain  Condition 
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Stress  (MPa) 


Figure  6.5.  Finite  Element  and  FIDEP  Predictions  for  the  Stresses  at  the  Interface  in 
Ti-24A1-1  INb  Matrix  After  Cool-Down  from  Processing  Temperature  of 
1010°C. 
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Plastic  Strain  (mm/mm) 


Figure  6.6  Finite  Element  and  FIDEP  Predictions  of  the  Plastic  Strains  in  Ti-24A1-1  INb 
Matrix  at  the  Fiber-Matrix  Interface  After  Cool-Down  from  Processing 
Temperature  of  1010°C. 
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Stress  (MPa) 


Plastic  Strain  (mm/mm) 


Figure  6.8  Variation  of  the  Plastic  Strains  Across  the  Cross-Section  of  the  CCM 
After  Cool-Down  from  Processing  Temperature  of  1010°C. 


the  plot,  are  compared  with  results  for  the  same  problem  solved  by  FEM  represented  by  the 
solid  lines.  The  agreement  with  the  finite  element  results  is  excellent.  The  computations 
show  that  the  matrix  material  reaches  the  yield  surface  at  approximately  600°C  and  remains 
in  contact  with  the  yield  surface  for  all  lower  temperatures.  The  axial  and  hoop  stresses  are 
tension  and  the  radial  stress  is  compressive  at  room  temperature  and  of  equal  magnitude. 

The  plastic  strains  in  the  matrix  at  the  fiber-matrix  interface  as  a  function  of 
temperature  are  shown  in  Fig.  6.6.  It  can  be  seen  here,  as  in  Fig.  6.5,  that  the  matrix 
deforms  elastically  until  600°C  where  the  matrix  stresses  reach  yield  surface  and  start 
deforming  plastically  after  which  measurable  plastic  strains  develop  in  the  three  principal 
directions.  These  results  are  compared  with  FEM  computations  and  provide  excellent 
agreement  as  shown  in  the  figure. 

Figures  6.7-6.8  display  the  stresses  and  the  plastic  strains  at  room  temperature 
across  the  cross-section  of  the  composite,  respectively.  The  shaded  areas  denote  the  regions 
of  plastic  deformation.  Plasticity  initiates  at  the  interface  and  propagates  across  the  matrix 
with  increased  mechanical  stress  on  the  matrix.  In  this  particular  case,  the  plastic 
deformation  has  stopped  at  r/b=0.91  and  the  rest  of  the  matrix  has  remained  elastic.  If  a 
tensile  mechanical  load  is  superimposed  at  room  temperature  this  will  further  propagate  the 
plastic  region  until  all  the  matrix  yields. 

6.3.2  Cyclic  Thermal  Loading 

In  this  example  the  SCS-6/Ti-24Al- 1  INb  composite  is  thermally  cycled  between 
processing  temperature  and  room  temperature  as  shown  in  Fig.  6.9.  The  effective  stress 
and  the  radial  stress  as  a  function  of  temperature  are  plotted  and  compared  with  results  from 
FEM  in  Fig.  6.10.  FfDEP  results  are  in  excellent  agreement  with  finite  element  results 
prior  to  attaining  a  temperature  of  850°C  at  the  end  of  the  first  cycle.  The  two  methods 
digress  after  850°C  because  of  the  changing  of  the  compressive  matrix  radial  stresses  at  the 
interface  to  tensile  at  this  temperature  (Fig.  6. 10).  In  the  case  of  the  tensile  radial  stress  the 
finite  element  model  asssumes  debonding  of  the  matrix  from  the  fiber  to  occur  whereas  the 
present  model  assumes  perfect  bonding.  Thus,  the  two  methods  predict  different  but  similar 
states  after  debonding  takes  place. 
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Loading  File  1 
Cyclic  Thermal  Loading 

Between  Room  Temperature  and  Processing  Temperature 
Temperature:  20  C  to  1010  C 
4 

0  30  90  120 

1010  20  1010  20 

0  0  0  0 


(a) 


Temperature  (°C) 


Computational  Steps 

(b) 


Figure  6.9  a)  Loading  Input  File,  and  b)  Schematic  of  the  Temperature  History  for 
Thermal  Cyclic  Loading  Simulation  of  SCS-6/Ti-24Al-l  INb 
Unidirectional  Composite. 
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Stress  (MPa) 


Figure  6.10  F1DEP  and  FEM  Predictions  for  the  Effective  and  Radial  Stresses 

in  the  Matrix  at  the  Fiber-Matrix  Interface  for  Thermal  Loading  History 
shown  in  Fig.  6.9. 
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6.4  Thermomechanical  Fatigue  (TMF) 

The  FIDEP  code  was  used  to  compute  stresses  for  simulated  TMF  cycles  to 
illustrate  the  micromechanical  stresses  developed  during  two  typical  TMF  tests  on  a  SCS- 
6/Ti-24Al-l  INb  composite.  The  two  tests  are  in-phase  (IP)  and  out-of-phase  (OOP)  cycles 
and  are  shown  schematically  in  Fig.  6.11.  In  both  tests,  the  first  part  of  the  cycle  is  the 
cool-down  from  processing  temperature  discussed  previously.  This  is  accomplished  under 
zero  external  load  conditions.  The  horizontal  axis  in  Fig.  6.1 1  represents  an  arbitrary  time 
scale  and  is  shown  in  terms  of  computational  steps.  Step  20  represents  the  room 
temperature  condition  after  cool-down.  The  values  shown  in  the  figure  correspond  to  TMF 
cycles  which  were  evaluated  experimentally  and  involve  a  maximum  stress  of  800  MPa, 
R=0. 1  (ratio  of  minimum  to  maximum  stress)  for  the  OOP  and  IP  cases,  and  a  temperature 
range  of  150  to  650°C. 

The  variation  of  the  stresses  across  the  cross  section  for  an  OOP  cycle  is  illustrated 
in  Fig.  6.12.  The  stresses  in  both  fiber  and  matrix  are  shown  at  the  minimum  temperature 
(maximum  load)  in  (a)  while  stresses  at  maximum  temperature  (minimum  load)  are 
presented  in  (b).  At  150°C,  Fig.  6.12(a)  shows  a  large  tensile  axial  stress  in  the  fiber  and 
tensile  axial  stresses  in  the  matrix.  The  fiber  stresses  are  essentially  independent  of  radial 
coordinate,  but  the  matrix  stresses  vary  significantly  with  radius  because  of  the  complex 
three-dimensional  stress  state  which  arises  in  the  concentric  cylinder  configuration.  Matrix 
axial  stresses  are  maximum  at  the  outer  radius  and  minimum  at  the  fiber-matrix  interface. 
Hoop  stresses,  on  the  other  hand,  are  minimum  at  the  outer  radius  of  the  cylinder  but  do  not 
vary  as  much  as  the  axial  component.  The  matrix  radial  stress  is  negative  at  the  fiber-matrix 
interface  and  goes  to  zero  at  the  outer  radius  because  of  the  imposed  stress-free  boundary 
condition.  The  stresses  developed  are  the  net  result  of  the  thermal  stresses  developed  at  a 
low  temperature,  the  axial  component  of  which  is  compression  in  the  fiber  and  tension  in  the 
matrix,  and  the  applied  mechanical  stresses  which  are  tension  in  both  fiber  and  matrix  in  the 
axial  direction. 

At  elevated  temperature  in  an  OOP  TMF  cycle,  the  stress  state  is  generally  small  as 
shown  in  Fig.  6.12(b).  Under  this  condition,  the  thermal  stresses  are  small  because  at  high 
temperature  the  stresses  are  relaxed  while  the  mechanical  stresses  are  small  because  this  is 
the  condition  of  minimum  load. 
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Cool-Down  and  In-Phase  TMF  Loading 
Load:  80  MPa  to  800  MPa 
Temperature:  150  C  to  650  C 
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Figure  6. 1 1  a)  Loading  Input  Files,  and  b)  Schematic  of  the  Temperature  and  Stress 
History  for  Typical  In-Phase  and  Out-of-Phase  TMF  Cycles. 
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Temperature  (°C) 


Stress  (MPa) 


Figure  6. 1 2a  Variation  of  Predicted  Stresses  with  Radius  for  an  Out-of-Phase  TMF 
Cycle  at  150°C. 


5i 


Stress  (MPa) 


1200 
1000 
800 
600 
400 
200 
0 

-200 
-400 

0  0.2  0.4  0.6  0.8  1 

Radius  (r/b) 

(b) 


- 

L  Fiber 

Matrix 

7  - Radial  Stress 

-  - Hoop  Stress 

;  — e —  Axial  Stress 

^  Fiber/Matrix 

- 

^  Interface 

1 

_ i _ i _ i _ i _ i _ i  t  i _ j _ i _ i _ 

i _ i _ i _ i _ t _ i _ i _ i _ 

Figure  6.12b  Variation  of  Predicted  Stresses  with  Radius  for  an  Out-of-Phase  TMF  Cycle 
at  650°C. 
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The  axial  stresses  in  the  matrix  at  the  fiber-matrix  interface  are  shown  in  Fig.  6.13. 
The  steps  on  the  horizontal  axis  correspond  to  those  of  the  TMF  spectrum  schematic  in  Fig. 
6.1  lb.  As  noted  in  the  discussion  of  the  previous  figure,  the  OOP  cycle  has  tensile  thermal 
stresses  as  well  as  tensile  mechanical  stresses  at  minimum  temperature  which  corresponds 
to  maximum  load.  This  results  in  a  large  stress  range  in  the  matrix  from  approximately  450 
MPa  at  minimum  temperature  (step  60)  to  approximately  -50  MPa  at  maximum  temperature 
(step  80).  The  IP  cycle,  on  the  other  hand,  has  a  much  smaller  stress  range  in  the  matrix 
material.  At  minimum  temperature  (step  60),  there  is  a  thermal  residual  stress  but  little 
contribution  from  applied  (minimum)  load.  The  net  stress  is  approximately  275  MPa. 
When  maximum  temperature  is  reached  (step  80),  the  thermal  stress  decreases  and  the 
mechanical  stress  increases.  The  net  result  is  a  stress  decrease  to  approximately  200  MPa. 
In  summary,  the  axial  stress  range  in  the  matrix  is  large  in  an  OOP  cycle  and  small  in  an  IP 
cycle. 


In  a  similar  fashion,  the  axial  stress  peaks  in  the  fiber  are  compared  for  an  IP  and 
OOP  TMF  cycle  in  Fig.  6.14.  In  the  fiber,  the  thermal  stresses  at  low  temperature  are 
compressive  because  of  the  low  value  of  a  compared  to  that  in  the  matrix.  For  an  IP  cycle, 
therefore,  the  axial  stresses  range  from  compressive,  due  to  residual  thermal  stress  at  low 
temperature  (step  60),  to  tension  due  to  the  applied  maximum  load  at  high  temperature  (step 
80).  In  the  OOP  cycle,  the  compressive  residual  stress  at  low  temperature  is  offset  by  the 
larger  mechanical  stress  due  to  maximum  load,  resulting  in  a  net  tensile  stress  of 
approximately  1000  MPa  at  minimum  temperature  (step  60).  Decreasing  the  load  while 
increasing  the  temperature  reduces  the  mechanical  stress  but  also  relaxes  the  compressive 
residual  stress.  Since  the  mechanical  stress  range  exceeds  the  thermal  stress  range  for  this 
set  of  conditions,  the  net  effect  is  a  reduction  in  stress  to  approximately  200  MPa  at 
maximum  temperature  (step  80).  Thus,  an  in-phase  cycle  produces  a  larger  stress  range  and 
higher  maximum  stress  in  the  fiber  than  an  out-of-phase  cycle  for  the  conditions  evaluated 
here.  This  is  in  contrast  to  the  stress  range  in  the  matrix  which  is  maximum  in  the  OOP 
cycle  and  smaller  in  the  IP  cycle. 

The  radial  and  hoop  stresses  in  the  matrix  at  the  fiber-matrix  interface  are  plotted  in 
Fig.  6. 15.  Since  the  applied  load  is  in  the  axial  direction,  its  contribution  to  hoop  and  radial 
stresses  is  minimal  and  arises  solely  due  to  Poisson's  ratio  and  3-D  plasticity  effects.  Thus, 
the  stresses  are  mostly  due  to  the  thermal  cycling,  and  the  IP  and  OOP  cycles  produce  a 
similar  stress  state.  At  minimum  temperature  (step  60),  the  hoop  stresses  in  the  matrix  are 
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Axial  Stress  (MPa) 


Computation  Step 

Figure  6.13  Axial  Stress  Predictions  in  the  Ti-24A1-1  INb  Matrix  at  the  Fiber-Matrix 
Interface  for  In-Phase  and  Out-of-Phase  TMF  Cyclic  Loading. 
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Axial  Stress  (MPa) 


Figure  6.14  Axial  Stress  Peaks  Predicted  in  the  SCS-6  Fiber  for  In-Phase  and  Out-of-Phase 
TMF  Cyclic  Loading. 
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tensile  whereas  the  radial  stresses  are  compressive.  At  the  maximum  temperature  (step  80), 
both  components  relax  to  nearly  zero. 
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PROGRAM  FIDEP3B 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

r 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 


THIS  PROGRAM  COMPUTES  ELASTIC -PLASTIC  STRESSES  IN  AXISYMMETRIC 
CONCENTRIC  CYLINDERS  SUBJECTED  TO  THERMOMECHANICAL  CYCLIC 
LOADING. 

INPUT  FILES: 

FDMAT  TEMPERATURE  DEPENDENT  MATERIAL  DATA 
FDLOAD  THERMOMECHANICAL  CYCLIC  LOADING  PROFILE 
OUTPUT  FILES: 

AT  INTERFACE: 

FDOUTS  STRESSES  AND  AXIAL  STRAIN  AT  EACH  STEP 
FDOUTE  STRAINS  AT  EACH  STEP 

FDSUM  STRESSES  AND  STRAINS  CORRESPONDING  TO  FDLOAD 
AT  CROSS-SECTION: 

FDUSP  STRESSES  AND  STRAINS  AT  THE  END  OF  THE  PROGRAM 
AT  THE  COMPOSITE  CROSSSECTION 

FINAL  REVISION:  13-NOVEMBER-1991 
PROGRAMMER:  DEMIRKAN  COKER  (255-1361) 

SUPERVISOR:  NOEL  E.  ASHBAUGH 

UNIVERSITY  OF  DAYTON  RESEARCH  INSTITUTE 


********************************************************************** 

DEFINITION  OF  THE  VARIABLES  USED  IN  THE  PROGRAM 
/MAT1/ 

E(201)  MODULUS  AT  RADIUS  R 

VMU(201)  POISSON’S  RATIO  AT  EACH  R 

CTE(201)  THERMAL  EXPANSION  COEFFICIENT  AT  EACH  R 

R(201)  RADIUS 

VF  FIBER  VOI  IJMF  FRACTION 

EZ  EZ  =  EZSTAR  -  2/EC*SR(INT+l)*VF*(VMU(INT)-VMU(INT+l)) 

I.E.  TOTAL  AXIAL  STRAIN  =  EZTOT 

/MAT2/ 

EC  RULE  OF  MIXTURES  COMPOSITE  MODULUS 
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C  TEMP(5,20) 

C 

C  EMAT(5 ,20) 

C  CTEMAT(5,20) 
C  VMUMAT(5,20) 
C  /MAT3/ 

C  DTC200) 

C  T(200) 

C  ET(5,200) 

C 

C  VMUT(5,200) 

C  CTET(5,200) 

C  PZC200) 

C  /MAT4/ 

C  TINI.TFIN 
C  PINI.PFIN 
C 

C  TREF 
C 

C  /MAT5/ 

C  SYMAT(5,20) 

C  EPLMAT(5,20) 
C 

C  /MAT6/ 

C  EPLT(5,200) 

C 

C  SY(5,200) 

C 

C  /LIMITS/ 

C  NTOT 
C  NRA 
C  INT 

C  /STRAINS/ 

C  ERTOT(201) 

C  ETTOT(201) 

C  EZTOTC201) 


TEMPERATURES  AT  WHICH  PROPERTIES  ARE  GIVEN  FOR 
CONSTITUENTS 

MODULI  AT  TEMPERATURE  FOR  THE  CONSTITUENTS 

CTES  AT  TEMPERATURE  FOR  THE  CONSTITUENTS 

POISSON'S  RATIOS  AT  TEMPERATURE  FOR  THE  CONSTITUENTS 

CURRENT  TEMPERATURE  MINUS  REFERENCE  TEMPERATURE 
CURRENT  TEMPERATURE 

CONSTITUENT  MODULI  INTERPOLATED  FOR  THE  CURRENT 
TEMPERATURE 

CONSTITUENT  POISSON'S  RATIIO  . 

CONSTITUENT  CTE  INTERPOLATED  . 

APPLIED  STRESS  CORRESPONDING  TO  THE  CURRENT  STEP 

INITIAL  AND  FINAL  TEMPS  READ  FROM  THE  LOADING  FILE 
INITIAL  AND  FINAL  APPLIED  STRESSES  READ  FROM  THE 
LOADING  FILE 

REFERENCE  TEMPERATURE  FOR  SECANT  CTE  =  PROCESSING 
TEMPERATURE  WHERE  THE  STRESSES  ARE  ASSUMED  TO  BE  ZERO 

ORIGINAL  YIELD  STRESS  AT  DEFINED  TEMPERATURES 
PLASTIC  MODULUS  FOR  THE  CONSTITUENTS  AT  DEFINED 
TEMPERATURES 

CONSTITUENT  PLASTIC  MODULI  INTERPOLATED  FOR  THE 
CURRENT  TEMP 

CONSTITUENT  ORIGINAL  YIELD  INTERPOLATED  FOR  THE 
CURRENT  TEMP 

NUMBER  OF  NODES  ALONG  THE  RADIUS  MINUS  1  =  NRA/2 . 

NUMBER  OF  ROWS  OF  MATRIX  A 

NODE  IN  THE  FIBER  AT  THE  INTERFACE 

TOTAL  RADIAL  STRAIN 
TOTAL  HOOP  STRAIN 
TOTAL  AXIAL  STRAIN 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


/XX/ 

ERP(201) , ETPC201) , EZP(201) 

OLD  PLASTIC  STRAINS  IN  THE  PLSTRAIN  .NE.  NEW  PLASTIC 
STRAINS  IN  SUBROUTINE  SOLVE  WHICH  DOES  NOT  USE  THIS 
COMMON  BLOCK 

ERP2C201),  ETP2(201),  EZP2(201) 

PLASTIC  STRAINS  USED  IN  SUBROUTINES  PLSTRAIN  AND  SOLVE 

/LEG/ 

JAR  TO  READ  NRA  FROM  SUBROUTINE  SOLVE  AT  THE  BEGINNING 

SRC201)  RADIAL  STRESS  ACROSS  THE  CROSSSECTION 

ST(201)  HOOP  STRESS  ACROSS  THE  CROSS-SECTION 

SZC201)  AXIAL  STRESS  ACROSS  THE  CROSS-SECTION 

NOMAT  NUMBER  OF  MATERIALS  -  FOR  NOW  2 

NROW(5)  NUMBER  OF  ROWS  OF  TEMPERATURE  DEPENDENT  DATA 

SE(201)  CURRENT  YIELD  SURFACE  -  USED  IN  SUBROUTINE  PLSTRAIN 

TE(20)  TEMPERATURE  READ  FROM  LOADING  FILE 

PAT(20)  APPLIED  STRESS  READ  FROM  LOADING  FILE 

NOT(20)  TOTAL  NUMBER  OF  STEPS  FOR  EACH  HALF  CYCLE 


*************************  ********************************************* 


COMMON/MAT0/E  F , EM , CTEF , CTEM 

COMMON/MATl/EC , E(201) ,VMU(201) , CTE(201) , R(201) , VF , EZ 
COMMON/MAT2/TEMP(5 , 20) , EMAT(5 , 20) , CTEMAT(5 , 20) , VMUMAT(5 , 20) 
COMMON/MAT3/DT(200) ,T(200) ,ET(5, 200) ,VMUT(5 , 200) , CTET(5 , 200) , PZ(200) 
COMMON/MAT4/TINI , TREF ,TFIN, PINI , PFIN 
COMMON/MAT5/SYMAT(5 , 20) ,  EPLMAT(5 , 20) 

COMMON/MAT6/EPLTC5 , 200) ,  SY(5 , 200) 

COMMON/L IMITS/NTOT ,  NRA,  INT 

COMMON/STRAINS/ERTOTC201),  ETTOT(201),  EZTOT(201) 

COMMON/LEG/JAR 

COMMON/EMECH/ETHERMAL , ERMECH , ETMECH , EZMECH 
REAL  ERP2(201) ,ETP2(201) ,EZP2(201) 

REAL  SR(201) ,  ST(201) ,  SZ(201),  SEFF(201),  SE(201) 

INTEGER  NOMAT,  NROW(5) 

INTEGER  NTIME 

REAL  TE(20) ,  PAT(20),  NOT(20) 


COMMON /WORK  S  P/RWK  S  P 
REAL  RWKSP( 100000) 

CALL  IWKINC100000) 

CALL  UMACHC-2,12) 

OPEN(UNIT=ll , FILE=  * FDMAT ' , STATUS= ’ OLD ' ) 
OPEN(UNIT=20 , FILE= ' FDLOAD ’ , STATUS= ’ OLD' ) 
OPEN(UNIT=12 ,  FILE='FDOUTS’ ,  STATUS=’NEW’) 
OPEN(UNIT=17,  FILE='FDOUTE ' ,  STATUS='NEW’) 
OPEN(UNIT=18 ,  FILE='FDUSP' ,  STATUS='NEW') 
OPEN(UNIT=19,  FILE=' FDSUM' ,  STATUS='NEW') 

NTIME  =  0 


C  -  -  -  Start  of  timed  sequence 
XT1=SECNDS(0.0) 


C  -  -  -  Read  dimension  NRA  of  the  matrix  A  to  be  solved 


JAR  =  0 

CALL  SOLVE(ILK , SR , ST, SZ , SEFF , B , ERP2 , ETP2 , EZP2) 
JAR  =  JAR  +  1 

PRINT  *,  'MATRIX  DIMENSION  IS',  NRA 
C  -  -  -  Read  loading  from  input  file  FDLOAD.DAT 


READ(20,*) 

READ(20,*)  ‘COMMENT  LINES  IN  INPUT  DATA 

READ(20,*) 

READ(20,*) 

READ(20,*)  NLOAD 

READ(20,*)  (NOT(I) ,  1=1, NLOAD) 

READ(20,*)  (TE(I),  1=1, NLOAD) 

READ(20,*)  (PAT(I) ,  1=1, NLOAD) 
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WRITE(12,*)  ’NO  OF  DATA  NLOAD 
WRITE(12,*)  'STEPS’,  (NOT(I),  1=1, NLOAD) 
WRITE(12,*)  'TEMP  ’,  (TE(I),  1=1, NLOAD) 
WRITEC12,*)  ’STRESS  ’,  CPAT(I),  1=1, NLOAD) 
TREE  =  TE(1) 


C  -  -  -  Read  temperature  dependent  material  properties 


CALL  RE ADMAT (NOMAT , NROW) 


C  -  -  -  Define  nodes  in  the  concentric  cylinder  model 
CALL  GEOMETRY (VF.B, A, DR, R) 


C  -  -  -  Print  header  lines 


CALL  HEADER(NTIME ,TE , PAT, NOT) 

CLOSE(ll) 

CLOSEC20) 


C  -  -  -  Start  cyclic  loading 


DO  995  ICYC  =  1,  NLOAD-1 

NSTEPS  =  N0T(ICYC+1)-N0T(TCYC) 
TINI  =  TE(ICYC) 

TFIN  =  TECICYC+1) 

PINI  =  PAT(ICYC)*1.0E6 
PFIN  =  PAT(ICYC+1)*1.0E6 

CALL  LOADING(NSTEPS, NOMAT) 

DO  990  ILK=2 ,NSTEPS+1 

NTIME  =  NTIME  +  1 
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C  -  -  -  Define  properties  at  each  node,  i.  e.  translate  properties  at 
C  temperature  for  material  1  and  2  into  node  properties;  E, alpha, nu 

C  In  addition,  define  shorthands  for  moduli,  Ef,  Em,  Ec. 

CALL  PROPERTIES (ILK) 

C  -  -  -  Calculate  stresses  with  old  plastic  strains  (dep-new  =  0) 

C  i.e.  assume  loading  increment  to  be  elastic  only 

CALL  SOLVE (ILK , SR , ST, SZ , SEFF , B , ERP2 , ETP2 , EZP2) 

C  -  -  -  Iterate  to  find  plastic  strain  increments 

CALL  PLSTRAIN(ILK,SR,ST,SZ,SEFF,B,ERP2,ETP2,EZP2,SE) 

C  -  -  -  Output  stresses  and  strains 

CALL  OUTS(SEFF , SR , ST , SZ , ILK , NTIME , ERP2 , ETP2 , EZP2) 

990  CONTINUE 

CALL  OUTSUM(INT, SZ , SR , ST, SEFF , NTIME , TFIN , PFIN) 

995  CONTINUE 

C  -  -  -  Output  stresses  and  strains  at  the  cross-section 

CALL  OUTUSP(SEFF , SR , ST , SZ , ERP2 , ETP2 , EZP2 , NTOT) 

C  -  -  -  Print  out  runtime  of  the  code 

CALL  REALTIME(NRA, NTIME, XT1) 

99  F0RMAT(2X,I4, 1X,F7. 1,1X,F8 . 1,1X,5(F8 . 1,1X) ,E10. 3) 

STOP 

END 


f>6 


c 

c 

c 


SUBROUTINES 


C  - 

C  READ  INPUT  MATERIAL  PROPERTIES  FOR  THE  FIBER  AND  THE  MATRIX 

C  - 


SUBROUTINE  READMAT (NOMAT , NROW) 

COMMON/MAT2/TEMP(5 , 20) , EMAT(5 , 20) , CTEMATO , 20) , VMUMAT(5 , 20) 
C0MM0N/MAT4/TINI ,TREF ,TFIN ,PINI ,PFIN 
C0MM0N/MAT5/SYMAT (5 , 20) ,  EPLMAT(5 , 20) 

REAL  TEMP,EMAT,CTEMAT,VMUMAT,TREF,TFIN>TREF2(5) 

INTEGER  NOMAT, NROWC5) 


C  -  -  -  Read  temperature  dependent  material  data 


READCll,*) 

READCll,*) 

READ(11,*)N0MAT 
DO  144  1=1, NOMAT 
READCll,*)  NROW(I) 

READCll,*) 

DO  144  J=l,NROW(I) 

READCll,*)  TEMP(I, J),EMAT(I , J),CTEMAT(I , J) .VMUMATCI , 3) , 

+  SYMAT(I , J) ,  EPLMATCI.J) 

WRITE(12 ,99)1 , J ,TEMP(I , J) ,EMATCI , J) ,CTEMAT(I , J) ,VMUMAT(I , J) , 
+  SYMAT(I.J),  EPLMAT(I ,3) 

144  CONTINUE 

99  FORMAT(2(I5) ,2X,4(F8 . 2) ,E9 . 2 ,F8 . 2) 

READCll,*) 

READCll, *)TREF 

RETURN 

END 


bl 


c  - 

C  COMPUTE  RADIUS  OF  NODES  AND  INTERFACE  NODE  IN  THE  CCM 

C  - 

SUBROUTINE  GEOMETRY(VF , B , A , DR , R) 

COMMON/L IMITS/NTOT ,  NRA,  INT 
REAL  VF,  B,  A,  DR,  R(201) 

INTEGER  NTOT,  INT 

READ(11,*) 

READ(11,*,END=44)  VF,  B 
WRITE(12,*)'VF=  VF 
GOTO  45 

44  VF  =  0.35 
B  =  1. 

45  A  =  SQRT(VF)/B 

PRINT  *,  'VF,  B  ' ,  VF,  B 
NTOT  =  NRA/2.0 
C  -  -  -  Nodes  in  Fiber 
R(l)  =  0.02E-10 
R(2)  =  A/2 
R(3)  =  A  -  A/100.0 
C  -  -  -  Nodes  in  Matrix 

R(4)  =  A  +  A/100.0 
DR  =  (B-A)/(NTOT-3) 

DO  10  1=5 ,NTOT+l 

R(I)  =  A  +  DR+CI-4) 

10  CONTINUE 

INT  =  3  !  A/DR  +  1 

DO  I=l,NTOT+l 
PRINT  *,  I,  R(I) 

END  DO 

PRINT  *,  ’  rf,  rm,  A  R(INT),  R(INT+1),A 
C  RCINT+1)  =  A 
RETURN 
END 
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c  - 

C  COMPUTE  INCREMENTAL  LOADING  AND  INTERMEDIATE  PROPERTIES 

C  - 


SUBROUTINE  LOADING(NSTEPS , NOMAT) 

COMMON/MAT3/DT (200) , T(200) , ET(5 , 200) , VMUT(5 , 200) , CTET(5 , 200) , PZC200) 
COMMON /MA T4/TI N I ,TREF ,TFIN , PINI , PFIN 

COMMON/MAT2/TEMP(5 , 20) , EMAT(5 , 20) , CTEMAT(5 , 20) , VMUMAT(5 , 20) 
COMMON/MAT5/SYMATC5 ,20) ,  EPLMAT(5 ,20) 

COMMON/MAT6/EPLT(5 , 200) ,  SY(5 , 200) 

DO  110  I=1,NSTEPS+1 

DTI I  =  (I-1.)*(TFIN-TREF)/NSTEPS/1. 

DTEMP  =  (I-1.)*(TFIN-TINI)/NSTEPS/1. 

T(I)  =  TINI  +  DTEMP 

DPZ  =  (I-1.)*(PFIN-PINI)/NSTEPS/1. 

PZ(I)  =  PINI  +  DPZ 
DT(I)  «  T(I)  -  TREF 
TT  =  TCI) 


C  -  -  -  Compute  the  material  properties  at  each  loading  temperature 


DO  K=l, NOMAT 
J=1 

DO  WHILE(TT. GT.TEMPCK, J)) 

J=J+1 
END  DO 

TFACTOR=(TT-TEMP(K, J-1))/(TEMP(K , J)-TEMP(K , J-l)) 
CALL  INTERP(I ,K,J ,TFACTOR,ET,EMAT) 

CALL  INTERPCI , K , J , TFACTOR , CTET , CTEMAT) 

CALL  INTERPCI, K,J, TFACTOR, VMUT.VMUMAT) 

CALL  INTERPCI , K , J , TFACTOR , SY , SYMAT) 

CALL  INTERPCI, K,J, TFACTOR, EPLT.EPLMAT) 

END  DO 
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110 


CONTINUE 


DO  120  K=l, NOMAT 
DO  120  I=1,NSTEPS+1 
SY(K,I)  =  SY(K,I)*1.E6 
EPLT(K,I)  =  EPLT(K,I)*1.E9 
120  CONTINUE 
RETURN 
END 


C  - 

C  INTERPOLATION  SCHEME  USED  IN  DETERMINING  INTERMEDIATE  PROPERTIES 

C  - 


SUBROUTINE  INTERP(I,K,J,TFACTOR,VALT,VALMAT) 

REAL  VALT(5 , 200) ,  VALMAT(5 ,20) ,  TFACTOR 
INTEGER  I,  K,  J 

VALT(K ,I)~VALMAT(K , J - 1)+(VALMAT(K , J)-VALMAT(K , J-l))*TFACTOR 

RETURN 

END 


C  - 

C  ASSIGN  PROPERTIES  TO  EACH  NODE 

C  - 


SUBROUTINE  PROPERTIES(ILK) 

COMMON/MAT0/EF , EM , CTEF , CTEM 

COMMON/MATl/EC , E(201) ,VMU(201) , CTE(201) , R(Z01) ,VF , EZ 
COMMON/MAT3/DTC200) ,T(200) , ET(5 , 200) ,VMUT(5 , 200) , CTET(5 , 200) , PZ(200) 
COMMON/MA T6/EPLT(5 , 200) ,  SY(5 , 200) 

COMMON/LIMITS/NTOT ,  NRA ,  INT 


EF  =  ET(1,ILK)*1.E9 
EM  =  ET(2 , ILK)*1. E9 
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1414E9 
! 113 . ZE9 


CTEF  =  CTET(l,ILK)*1.0E-6 
CTEM  =  CTET(2 ,ILK)*1.0E-6 
DO  21  1=1, INT 

E(I)  =  ET(1,ILK)*1. E9 
VMU(I)  =  VMUT(l.ILK) 

CTECI)  =  CTET(l,ILK)*l.E-6 

21  CONTINUE 

DO  22  I=INT+l,NTOT+l 
E(I)  =  ET(2,ILK)*1.E9 
VMU(I)  =  VMUTC2.ILK) 

CTE(I)  =  CTET(2,ILK)*l.E-6 

22  CONTINUE 

EC  =  EF*VF  +  EM*(1-VF) 

RETURN 

END 


I414.E9 
!0. 3 
14.7E-6 


1113.7E9 

*0.22 

I9.44E-6 
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c  - 

C  COMPUTE  STRESSES  GIVEN  A  PLASTIC  STRAIN  INCREMENT 
C  - 


SUBROUTINE  SOLVE(ILK,SR,ST,SZ,SEFF,B,ERP,ETP,EZP) 

PARAMETER(NRA=100,  LDA=NRA ,  IPATH=1) 

COMMON/LIMITS/NTOT.NNRA ,  INT 
COMMON/MAT0/E F , EM , CTEF , CTEM 

COMMON/MAT1/EC , E(201) ,VMU(201) , CTE(201) , R(201) , VF ,  EZ 
COMMON/MAT3/DT (200) ,T(200) , ET(5 , 200) ,VMUT(5 , 200) , CTET(5 , 200) , PZ(200) 
COMMON/MAT6/EPLT(5 , 200) ,  SY(5 , 200) 

REAL  ERP(201),  ETP(201),  EZP(201) 

COMMON/STRAINS/ERTOT(201),  ETTOT(201),  EZTOT(201) 

COMMON/LEG/JAR 

REAL  AA(201),  BBC201),  CC(201),  DD(201),  FFC201),  GGC201),  HHC201) 
REAL  QQC201),  PLAS(201),  PP(201) 

REAL  AMAT(NRA,NRA),  BMAT(NRA),  XSOL(NRA) 

REAL  SR(201),  ST(201),  SZ(201),  SEFF(201) 

NNRA  =  NRA 
IF  (JAR. EQ.0)THEN 
RETURN 
END  IF 


C  -  -  -  Calculate  Ez  without  the  radial  stress  component  at  the  interface 


TEGRAL  =  0. 

DO  I=IN1 +1 , NTOT+1 

TEGRAL=TEGRAL+(RCI)-RCI-l))*(EZP(I)*R(I)+EZP(I-l)*R(I-l))/2 
END  DO 

EZSTAR  =  PZ(ILK)/EC  +  DT(ILK)* 

(CTEF*EF*VF+CTEM*EM*(1-VF))/EC 
EZSTAR  =  EZSTAR  +  2*TEGRAL*EM/B**2/EC 


C  -  -  -  Compute  coefficients  for  the  FD  equations  to  use  in  the  A-matrix 
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DO  30  1=2 ,NTOT+l 

AA(I)  =  (R(I)-R(I-l))/2/R(I) 

BB(I)  =  E(I)/E(I-1) 

CC(I)  =  R(I)/R(I-1) 

DD(I)  =  (1+VMU(I))*(VMU(I)+AA(I)) 

FF(I)  =  C1+VMU(I))*(1-VMU(I)+AACI)) 

GG(I)  =  (1+VMU(I-1))*(VMU(I-1)-AA(I)*CC(I))*BBCI) 

HH(I)  =  (1+VMU(I-1))*(1-VMU(I-1)-AA(I)*CC(I))*BB(I) 

QQCI)  =  ECI)*DTCILK)*(CTECI)*C1+VMU(I))-CTE(I-1)*C1+VMUCI-1))) 
PLAS(I)  =  AA(I)*CERPCI)+CCCI)*ERPCI-1))  -  ETP(I)*(1+AA(I)) 

+  +  ETP(I-1)*(1-AA(I)*CC(I))  +  VMU(I)*(ERP(I)+ETP(I)) 

+  -  VMU(I-1)*(ERP(I-1)+ETP(I-1)) 

PLAS(I)  =  E(I)*PLAS(I) 

PP(I)  =  PLAS(I)  +  E(I)*EZSTAR*CVMUCI)-VMU(I-1)) 

30  CONTINUE 


C  -  -  -  Initialize  A  and  B 

DO  51  I=l,2*NTOT 
DO  52  J=l,2*NTOT 
AMAT(I,J)=  0.0 
52  CONTINUE 

BMAT(I)  =  0.0 
51  CONTINUE 

C  -  -  -  Construct  the  A-matrix  and  B-matrix 


AMAT(1, 1)  =  AAC2) 

AMAT(NTOT+l, 1)  =  (1  +  VMU(1))*(1  -  2*VMU(1))*BB(2) 
DO  55  1=2 ,NTOT 

AMAT(I-1,I)  =  -1.0  !  upper  left 

AMAT(I,I)  =  1/CCCI+l) 

AMAT(NT0T+I-1,I)  =  DD(I)  !  lower  left 

AMATCNTOT+1,1)  =  -GG(I+1) 

AMAT(I,NT0T+I-1)  =  AA(I+1)  !  upper  right 

AMAT(NTOT+I .NTOT+I-l)  =  HH(I+1)  !  lower  right 
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55  CONTINUE 


C  -  -  -  This  is  the  vmu  correction  at  the  interface  from  the  axial  strain 


AMAT(NTOT+INT,INT)  =  2*EM/EC*VF*(VMU(INT)-VMU(INT+1))**2  -  GG(INT+1) 
DO  66  1=1 , NTOT 

AMAT(I,NTOT+I)  =  AA(I+1)  !  uUyc,  right 

AMAT(NTOT+I , NTOT+I)  =  -FF(I+1)  !  lower  right 

BMAT(I)  =  0. 

BMAT(NTOT+I)  =  QQ(I+1)  -  PP(I+1) 

66  CONTINUE 


C  -  -  -  Solve  for  the  radial  and  hoop  stresses 


CALL  LSARG(NRA , AMAT , LDA , BMAT , IPATH ,XSOL) 


C  -  -  -  Enforce  boundary  conditions 


SRCNTOT+1)  =0.0 

ST(1)  =  XSOL(l)  !  XSOL(l)  =  SR(1) 


(  -  -  -  Extract  radial  and  hoop  stresses  from  XSOL 


DO  86  I=1,NT01 
SR(I)  =  XSOL(I) 

ST( I + 1)  =  XS0L(N10T+I) 

86  CONTINUE 

[/  -  E/STAR  2/F  *SR(INT+1)*VE *CVMU(INT)-VMU(INT+1)) 


Compute  axiat  and  effective  stresses 


\/)  i01  I  l.NTOWJ 

6/(1)  1 (I )♦(!/! /P( I ))  -  CTE(I)*E(I)*DT(IIK) 

.  VMH( I  )*( SRC  I  ).S1(I )) 

'  \ 


SQUARE  =  (SR(I)-ST(I))**2+(SR(I)-SZ(I))**2+(ST(I)-SZ(I))**2 
SEFF(I)  =  SQRT((SQUARE)/2.) 

101  CONTINUE 

C  -  -  -  Compute  total  strains 
DO  I=1,NT0T+1 

ERTOT(I)=(SR(I)-VMU(I)*(STCI)+SZCI)))/ECI>fCTE(I)*DT(ILK)+ERP(I) 

ETTOTCI)=CSTCI)-VMU(I)*CSR(I)+SZ(I)))/ECI)+CTECI)*DT(ILK)+ETP(I) 

EZTOT(I)=(SZ(I)-VMU(I)*(ST(I)+SR(I)))/E(I)+CTE(I)*DT(ILK)+EZP(I) 

END  DO 

RETURN 

END 

C  - 

C  SUBROUTINE  TO  COMPUTE  PLASTIC  STRAIN  INCREMENTS  USING  MENDELSON'S 
C  ALGORITHM 

C  - - - 

bUBROUTINE  PLSTRAIN(ILK,SR,ST,SZ,SEFF,B>ERP2,ETP2,EZP2,SE) 
COMMON/LIMITS/NTOT.NNRA,  INT 
COMMON/MAT0/EF , EM , CTEF , CTEM 

COMMON/MAT 1/EC , E(201) ,VMU(201) , CTE(201) , R(201) , VF , EZ 
COMMON/MAT3/DTC200) , T(200) , ET(5 , 200) ,VMUT(5 , 200) , CTET(5 , 200) , PZ(200) 
COMMON/MAT6/EPLT(5 , 200) ,  SY(5 , 200) 

COMMON/LEG/JAR 

COMMON/STRAINS/ERTOT(201) ,  ETTOT(201),  EZTOT(201) 

COMMON/XX/ERP(201) ,  ETP(201),  EZP(201) 

COMMON/XXFEB6/DRTEST1C201),  DRTEST2(Z01) ,  DRTEST3(201) 
COMMON/XXFEB6/DTTEST1C201),  DTTEST2(201) ,  DTTEST3(201) , EPEFF(201) 
REAL  SR(201),  ST(201),  SZ(201),  SEFF(201) 

REAL  SE(201) 

REAL  DERP(201),  DETP(201),  DEZP(201),  DEP(201) 

REAL  ERP2(201) ,  FTP2(201),  EZP2(201) 

REAL  ERPCHECK(201) ,  FTPCHECK(201) 


REAL  YLDFLGC201) 

REAL  DELTAR1(201),DELTAR2(201) 

REAL  DELTAT1(201) ,DELTAT2(201) 

REAL  ARCOEF(201),BRCOEF(201),ATCOEF(201),BTCOEF(201) 
INTEGER  K,  J,  JJJ 


C  DT :  Current  temperature  minus  the  reference  temperature 

ITER  =  0 
IZMIR  =  0 

C  -  -  -  Calculate  new  yield  surface  after  strain-hardening 

DO  I  *  INT+1,  NTOT+1 

XM  =  EPLT(2,ILK)/E(I) 

DSDE  =  XM*E(I)/(1.-XM) 

SE(I)  =  DSDE*EPEFF(I)  r  SY(2,ILK) 

END  DO 

C  -  -  -  Check  if  the  new  stresses  are  greater  than  new  Y.S. 

8  DO  I  =  INT+1,  NTOT+1 

IF(SEFF(I).GT.SE(I))  GO  TO  9 
END  DO 
RETURN 

C  -  -  -  If  there  is  yielding  of  some  nodes,  start  iterating  for  dep 

9  DO  55  KZ=1, 12 

DO  10  I  =  INT+1,  NTOT+1 

C  -  -  -  Calculate  new  yield  surface  due  to  plas  strain  incr.  at  node  i 

XM  =  EPLT(2 , ILK)/E(I) 

DSDE  =  XM*E(I)/(1.-XM) 

SECT)  =  DSDE*EPEFF(I)  +  SY(2,ILK) 
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F  =  SEFF(I)  -  SE(I) 


C  -  -  -  Check  stress  state  for  yielding  at  node  i 

IF(F.GT.0)  GOTO  50 
YLDFLG(I)  =  0. 

GOTO  10 

50  YLDFLG(I)  =2.0 

C  -  -  -  Calculate  modified  total  strains 

ERMTS  =  ERTOT(I)  -  ERP(I)  !  erp  is  w/o  derp 

ETMTS  =  ETTOT(I)  -  ETP(I)  !  W/O  DETP(I) 

EZP(I)  =  -ERP(I)  -  ETP(I) 

EZMTS  =  EZTOT(I)  -  EZP(I)  !  W/O  DEZP(I) 

551  «  ERMTS  -  ETMTS 

552  =  EZMTS  -  ETMTS 

553  =  ERMTS  -  EZMTS 

C  -  -  -  Calculate  equivalent  (effective)  modified  total  strain  (MTS) 

EET  =  SQRT(SS1**2+SSZ**2+SS3**2)*SQRT(2 . )/3 . 

C  -  -  -  Relationship  between  dep  or  PSI  and  eff  MTS 

DENOMDEP  =  1.  +  2./3.*(l+VMU(I))/E(I)*DSDE 
DEP(I)  =  EET  -  2./3.*(l+VMU(I))/E(I)*SE(I) 

DEP(I)  =  DEP(I)/DENOMDEP 

C  -  -  -  Calculate  new  PSIs  using  modified  Prandtl -Reuss  equations 

DERP(I)  =  DEP(I)/3. /EET* (2* ERMTS -ETMTS -EZMTS) 

DETP(I)  =  DEP(I)/3 ./EET* (2* ETMTS -ERMTS -EZMTS) 

DEZ2  =  DEP(I)/3 . /EET* (2* EZMTS -ETMTS -ERMTS) 

DEZP(I)  =  -(DERP(I)  +  DETP(I))  !  Compare  dezp  &  dez2 
ERP2(I)  =  ERP(I)  +  DERP(I) 
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ETP2(I)  =  ETP(I)  +  DETP(I) 

EZP2(I)  =  -ERP2(I)  -  ETP2CI) 

EZP2TEST  =  EZP(I)  +  DEZP(I) 

10  CONTINUE 

C  -  -  -  Solve  O.D.E.s  with  new  plastic  strain  increments 

11  CALL  SOLVE(ILK,SR,ST,SZ,SEFF,B,ERP2,ETP2,EZP2) 


C  -  -  -  Save  old  plastic  strains  from  the  three  previous  iterations 

DO  I  =  INT+1,  NTOT+1 
DRTESTl(I)  =  DRTEST2(I) 

DRTEST2(I)  =  DRTEST3CI) 

DRTEST3(I)  =  ERP2CI) 

DTTESTl(I)  =  DTTEST2CI) 

DTTEST2CI)  =  DTTEST3(I) 

DTTEST3CI)  =  ETP2CI) 

END  DO 


C  -  -  -  Faster  convergence  scheme  extrapolates  the  difference  between 
C  previous  plastic  strain  increments  to  zero. 


IZMIR  =  IZMIR  +  1 
IF (IZMIR . EQ . 4)THEN 
IZMIR  =  0 

DO  22  I  =  INT+1,  NTOT+1 
PRODTEST  =  DRTEST1(I)*DRTEST2(I)*DRTEST3(I) 

DRDRTEST  =  ABS(DRTEST2(I)-DRTEST3(I)) 

DTDTTEST  =  ABS(DTTEST2(I)-DTTEST3(I)) 

I F (PRODTEST . NE . 0 . 000000 . AND . DRDRTEST . GT . IE -8 . AND . 

+  DTDTTEST. GT. l.E-10)THEN 

DELTARl(I)  =  DRTEST2(I)  -  DRTESTl(I) 

DELTAR2(I)  =  DRTEST3(I)  -  DRTEST2(I) 

ARCOEF(I)  =  (DELTAR1(I)-DELTAR2(I))/(DRTEST2(I)-DRTEST3(I)) 
BRCOEF(I)  =  DELTAR2(I)  -ARCOEF(I)*DRTEST3(I) 
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if (atcoef (i) . gt . le-14)then 

ERP2(I)  =  -  BRCOEF(I)/ARCOEF(I) 
end  if 

DELTAT1CI)  =  DTTEST2CI)  -  DTTESTl(I) 

0ELTAT2CI)  =  DTTESTBCI)  -  DTTEST2(I) 

ATCOEF(I)  =  (DELTAT1(I)-DELTAT2(I))/(DTTEST2(I)-DTTEST3(I)) 
BTCOEF(I)  =  DELTAT2CI)  -ATCOEF(I)*DTTEST3(I) 
i f (atcoef (i ) . gt . le -14)then 

ETP2(I)  =  -  BTCOEF(I)/ATCOEF(I) 
end  if 

EZP2(I)  =  -  ERP2(I)  -  ETP2CI) 

END  IF 

22  CONTINUE 

GOTO  11 
END  IF 

55  CONTINUE 

DO  M  =  INT+1,  NTOT+1 
ERP(M)  =  ERP2(M) 

ETP(M)  =  ETP2(M) 

EZP(M)  =  -ERP(M)  -ETP(M) 

EPEFF(M)  =  EPEFF(M)  +  DEP(M) 

END  DO 
RETURN 
END 


79 


c  - 

C  HEADER  LINES  FOR  THE  OUTPUT 

C  - 


SUBROUTINE  HEADER(NTIME , TE , PAT , NOT) 

REAL  TE(20) ,  PAT(20),  NOT(20) 

WRITE(12 , 750) 

WRITE(19,751) 

WRITE(17,752) 

WRITE(18 ,753) 

750  F0RMAT(1X, 'STEP  ’/TEMP  '  ,'APPSTRESS  ','SEFF  '/  SR  ’/  ST 

S  '  SZ  '/  SY(MPA)  ' ,  'EZMECH') 

751  F0RMAT(1X,  'STEP  '/TEMP  ','APPSTRESS  ',’SZF  ','  SZM  '/  SEFFM  ', 

$  '  EEFFM  ’,'  EZMECH  ') 

752  F0RMAT(1X, '  STEP  '/TEMP  '/ERPM  '/ETPM  ','  EZPM  ’/  ERMECH  ', 

$  '  ETMECH  ','  EZMECH  ',  'ETHERMAL') 

753  F0RMAT(1X,’R/B  '  ,'SEFF  '/  SR  '/  ST  ', 

S  '  SZ  '  /  ERP  ' ,  'ETP  ’,  'EZP') 

RETURN 

END 


C  - - - 

C  PRINT  STRESS  OUTPUT  IN  FDOUTS.DAT  AND  STRAIN  OUTPUT  IN  FDOUTE.DAT 

C  - 


SUBROUTINE  OUTS(SEFF,SR,ST,SZ,ILK,NTIME,ERPZ,ETP2,EZP2) 
COMMON/MAT0/EF , EM , CTEF , CTEM 

COMMON/MAT1/EC ,E(201) ,VMU(201) ,CTE(201) ,R(201) ,VF ,EZ 
COMMON/MAT3/DT(200) , T(200) , ET(5 , 200) , VMUT(5 , 200) , CTET(5 , 200) , PZ(200) 

COMMON/MAT6/EPLTC5 , 200) ,  SY(5 , 200) 

COMMON/LIMI TS/NTOT,  NRA ,  INT 

COMMON/STRAINS/ERTOT(201),  ETTOT(201),  EZTOT(201) 
COMMON/EMECH/ETHERMAL, ERMECH, ETMECH, EZMECH 
REAL  ERP2(201),ETP2(201),EZP2(201) 
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REAL  SR(201),  ST(201),  SZ(201),  SEFF(201) 


ICVC  =  ICYC  +  1 
IF(ICYC.EQ.1)THEN 

WRITE(12 ,99)  NTIME ,  TE(1),  PAT(1)/1. E6,  SEFF(1)/1.F6, 

+  SR(1)/1.E6,  ST(1)/1.E6,  SZ(1)/1.E6,  SY(2,1)/1.E6, EZMECH 

END  IF 

551  =  SEFF(INT+1)/1.E6 

552  =  SR(INT+1)/1.E6 

553  =  ST(INT+1)/1.E6 

554  =  SZ(INT+1)/1.E6 

555  =  SY(2,ILK)/1.E6 
ETHERMAL  =  CTE(INT+1)*DT(ILK) 

ERMECH  =  ERT0T(INT+1)  -  ETHERMAL 
ETMECH  =  ETTOTCINT+1)  -  ETHERMAL 
EZMECH  =  EZTOTCINT+1)  -  ETHERMAL 

WRITE(12,99)  NTIME,  T(ILK),  PZ(ILK)/1. E6,  SSI,  SS2,  SS3,  SS4, 
+  SS5,  EZMECH 
L  =  INT+1 

WRITE(17,101)  NTIME,  T(ILK),  ERPZ(L),  ETP2(L),  EZP2(L) 

+  ,  ERMECH,  ETMECH,  EZMECH,  ETHERMAL 

99  F0RMAT(2X,I4, 1X,F7. 1, 1X,F8 . 1,1X,5(F8 . 1, IX) ,E10. 3) 

101  F0RMAT(2X , 15 , 2X , F7 . 2 , 7(2X , Ell . 5)) 

RETURN 

END 
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c  - 

C  PRINT  SUMMARIZED  OUTPUT  IN  FDSUM.DAT 
C  - 


SUBROUTINE  OUTSUM(INT,SZ,SR,ST,SEFF,NTIME,TFIN,PFIN) 
COMMON/STRAINS/ERTOT(201),  ETTOT(201),  EZTOT(201) 
COMMON/EMECH/ETHERMAL , ERMECH , ETMECH , EZMECH 
REAL  SR(201),  ST(201),  SZ(201),  SEFF(201) 

XSZF  =  SZ(INT)/1.E6 
XSZM  =  SZ(INT+1)/1.E6 
XSEFFM  =  SEFrCINT+l)/l.E6 
IZ  =  INT+1 

SQUARE  =  (ERTOT(IZ)-ETTOT(IZ))**2  +  (ERTOT(IZ)-EZTOT(IZ))**2  + 
+  (ETTOT(IZ)-EZTOT(IZ))**2 

XEEFFM  =  SQRT(( SQUARE) *2 . j/3 . 

XEMECH  =  EZTOT(INT+l)  -  ETHERMAL 
WRITE(19,222)  NTIME ,  TFIN,  PFIN/1.E6, 

+  XSZF,  XSZM,  XSEFFM,  XEEFFM,  XEMECH 

222  FORMAT(I6 ,2(1X,F9.2),3(1X,F9.2),2(1X,E12.4)) 

RETURN 

END 


C  - 

C  PRINT  OUTPUT  AT  THE  CROSS-SECTION  IN  FDUSP.DAT 
C  - 


SUBROUTINE  OUTUSP(SEFF , SR , ST, SZ , ERP2 , ETP2 , EZP2 , NTOT) 
COMMON/MAT0/EF , EM , CTE  F , CTEM 

COMMON/MAT1/EC , E(201) , VMU(201) , CTE(201) , R(201) ,VF , EZ 
REAL  ERP2(201),ETP2(201),FZP2(Z01) 

REAL  SR(201),  ST(201),  SZ(201),  SEFF(201) 


DO  I=l,NTOT+l 

551  =  SEFF(I)/1.E6 

552  =  SRCI)/1.E6 
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553  =  ST(I)/1.E6 

554  =  SZ(I)/1.E6 

WRITE(18 ,88)  R(I) , SSI, SS2 ,SS3 , SS4,ERP2(I) , ETP2(I) , EZP2(I) 
END  DO 

88  FORMAT(2X , F7 . 4 , 2X , 4(F9 . 2 , IX) , 3(F9 . 5 , IX)) 

RETURN 

END 

C  - 

C  COMPUTE  RUNTIME 

C  - 


SUBROUTINE  REALTIME(NRA , NTIME , XT1) 

DELTA  -  SECNDS(XTl)/60. 

WRITE(12,*)  'RUNTIME  IS  DELTA,'  MIN  FOR  ',  NRA , '  NODES  ’ 
$  NTIME,  'TIME  INCRS.’ 

PRINT  *,  'RUNTIME  IS  '.DELTA,'  MIN  FOR  ',  NRA,'  NODES  ', 

$  NTIME,  'TIME  INCRS.' 

RETURN 

END 
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APPENDIX  B 

ELASTIC  SOLUTION  OF  TWO  CONCENTRIC  CYLINDERS 


Stresses  in  the  matrix  and  the  fiber  are  given  by  [28]; 


<Jrm  =  A(l-~)  +  Pr~, 
r  r 


r  r 


<x  =  C, 

zm  ’ 
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APPENDIX  C 
SAMPLE  OUTPUT  FILES 


8.5 


FDOUTS.DAT 
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